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1985-1986  AIR  FORCE  GEOPHYSICS  SCHOLAR  PROGRAM 


INTRODUCTION 

The  Geophysics  Scholar  Program  was  initiated  to  provide  research 
scholars  with  one  year  appointments  to  conduct  research  at  the  Air  Force 
Geophysics  Laboratory,  Hanscom  AFB,  MA. 

Extensive  mailings  were  made  to  technical  departments  at  universi¬ 
ties  around  the  United  States  where  programs  of  prime  interest  to  the 

Geophysics  Laboratory  were  established.  These  included  atmospheric  stud¬ 
ies,  space  science,  geophysics,  meteorology  and  related  applied 
sciences. 

Nine  scholars  were  appointed  during  the  period  starting  September 
1985  and  running  through  the  end  of  the  Contract  on  31  August  1986.  Two 
of  these  were  continued  in  this  Geophysics  Scholar  Program  for  a  second 
year  after  having  held  scholar  appointments  under  a  previous  contract. 

The  scholars  attended  7  meetings  or  conferences  during  their 
appointment  period.  Several  technical  papers  were  presented  by  the 

scholars  during  the  year.  The  final  technical  reports  on  the  scholar's 
work  are  included  in  this  report. 

This  program  was  judged  a  success  by  both  the  scholars  and  their 
laboratory  associates.  The  opportunity  of  having  new  research  people  on 
a  short  term  basis  was  felt  to  be  very  stimulating  and  worth  while. 
Their  interactions  with  the  laboratory  were  very  positive. 

The  program  was  scheduled  to  end  during  the  summer  of  1985  but  was 
extended  to  31  August  1986  in  order  to  accommodate  scholar  continuation 

of  2nd  year  selectees  from  the  1984-85  program.  As  a  result,  there  were 

only  9  scholars  during  the  1985-86  academic  year. 
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COLLISION  INDUCED  ELECTRONIC  TRANSITIONS  IN  N?+ 

by 

Anthony  V.  Dentamaro 

ABSTRACT 

An  optical -optical  double  resonance  technique  is  used  to  determine 
propensities  for  collision  induced  electronic  relaxation  by  helium  atoms 
from  a  specific  A^nu^(v=4,  J)  rotational  level  to  the  X^^g+(v=7)  manifold 
of  N2+.  The  propensities  for  colli sional  transfer  from  this  specific  level 
to  the  nearly  degenerate  (  ~  0.04  cm"l  separation)  spin  components  of  the 
X(v= 7)  state  are  resolved  by  scanning  the  probe  laser  through  the  B2eu+  - 
X^g4-  (5,7)  band  whose  upper  state  Is  perturbed.  The  results  show  the 
propensities  to  be  quite  different  and  strongly  dependent  on  the  A(v=4,J) 
level  initially  populated  by  the  pump  lasers.  The  observation  of  these 
propensities  for  collisional  electronic  energy  transfer  through  a  large  energy 
gap  of  approximately  1760  cm-*  demonstrates  the  remarkable  fact  that  this 
process  occurs  as  fast  or  faster  than  rotational  energy  transfer  through 
gaps  of  ~  10  cm~l.  These  results  are  found  to  be  in  qualitative  agreement 
with  theoretical  relative  cross-sections  derived  by  Alexander  and  Corey 
[J.  Chem.  Phys.  84,  100(1986)]  for  inelastic  collision  induced  transitions 
between  2n  and  2S  electronic  states  of  a  diatomic  molecule. 


I.  Introduction 


Selection  rules  for  optical  transitions  in  diatomic  molecules  are  well 
documented  and  listed  in  the  standard  texts  of  molecular  physics.  However,  very 
little  is  known  even  now  about  transitions  in  molecules  due  to  collisions  with 
atoms  or  other  molecules.  Although  such  collisions  have  been  examined  both 
theoretically  and  experimentally  for  many  years,  only  recently  have  there  been 
comprehensive  studies  of  the  resulting  propensity  rules.  The  majority  of  these 
experiments  and  models  have  been  concerned  with  vibrational  and  rotational 
energy  transfer,  and  very  little  has  been  done  in  the  line  of  collision- 
induced  electronic  transitions. 

Alexander  and  Corey2  have  expanded  upon  past  work  involving  rotational 
transition  propensity  rules  and  have  developed  a  quantum  mechanical  theory 
for  collisional  transitions.  Predictions  from  earlier  papers 

authored  by  Alexander3*4*5  have  been  experimentally  tested  with  the  results  being 
in  generally  good  agreement  with  the  theory.  In  particular,  Daniel  H.  Katayama 
of  AFGL  has  performed  many  of  these  tests,  and  during  this  past  year  we  have 
examined  the  propensities  for  transitions  in  the  A-X  Meinel  system  of  N2+ 
when  the  collision  partner  is  atomic  helium5. 

Both  the  quantum  formulation  and  the  experimental  results  show  under¬ 
lying  dipole  behavior  in  that  there  is  a  propensity  toward  aJ=0,±1  in 
these  electronic  transitions.  Aside  from  this  similarity,  the  agreement  between 
theory  and  experiment  is  only  fair.  This  is  perhaps  due  to  the  fact  that 
the  quantum  mechanical  formulation  does  not  explicitly  depend  on  the  form  of 
the  atom-molecule  potential.  Therefore,  the  next  step  would  be  to  calculate 
cross  sections  for  these  transitions  using  a  specific  interaction  between 
the  N2+  and  He.  Several  possible  forms  exist  for  the  missing  potential.  The 


dipole  behavior  suggests  a  long-range  attractive  interaction,  while  the 
deviation  from  strict  optical  selection  rules  may  indicate  the  existence 
of  repulsive  forces.  Of  particular  interest  is  the  collision  complex  model 
of  bound  atom-molecule  states  which  allows  for  the  inclusion  of  the  energy 
gap.  Thus  far,  collision  complex  models  have  been  considered  mainly  for 
vibrational  predissociation,  but  we  feel  that  an  electronic  analogue  is 
possible. 

Cross  sections  and  transition  rates  may  be  calculated  numerically 
for  any  of  these  models,  and  our  goal  is  to  find  the  right  potential  to 
reproduce  our  recent  experimental  results. 


1 1 .  Objectives  of  the  Research  Effort 

The  primary  goal  of  this  past  research  period  was  to  study  collision- 
induced  transitions  in  diatomic  molecules  and  determine  the  resulting 
propensities.  Specifically,  we  were  interested  in  electronic  transitions. 
Continuing  wortc  by  Daniel  H.  Katayama  at  AFGL  involved  using  an  optical  - 
optical  double  resonance  technique  in  order  to  obtain  rotational 
resolution  in  collision  experiments 

Various  theories  of  atom-molecule  collisions  exist  which  describe  the 
different  aspects  of  these  interactions  and  yield  a  wide  range  of  predictions 
for  combinations  of  molecules,  partners  and  energy  levels;  however,  there  has 
been  no  single  general  model  which  describe  all  of  these  systems.  Certainly, 
no  one  approach  is  able  to  account  for  our  observations.  Incorporating  the 
ideas  of  all  of  these  models,  our  next  objective  was  to  choose  a  specific 
potential  which  was  representati ve  of  the  atom-molecule  interaction  and 
calculate  the  resulting  cross  sections  for  comparison  with  our  experimental 
results  for  the  N^-He  system. 
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III.  Experimental  Method  and  Results 

The  two  laser,  double  resonance  system  used  to  conduct  the  experiment  is 
essentially  the  same  as  that  described  previously?*®.  The  "pump"  and  "probe" 
dye  lasers  are  excited  simultaneously  by  the  green  and  UV  light  pulses,  re¬ 
spectively,  of  a  Nd:YAG  laser.  The  pump  laser  with  Rh  640  dye  solution  se¬ 
lectively  populates  a  specific  rotational  level  of  the  N2  +  A2nin--X27g+ 

(4,0)  band.  The  collision  induced  electronic  transitions  by  helium  atoms  to  the 
X2Eg+(v=7)  rotational  manifold  which  lies  approximately  1760  cm-1  lower 
in  energy  relative  to  the  A(v=4)  level  are  detected  by  scanning  the  probe 
laser  with  Coumarin  450  dye  solution  through  the  B  *ru+-X2z;g+  (5,7) 
band.  The  N2+  are  formed  by  interacting  nitrogen  with  helium  metastable 
atoms  downstream  from  a  dc  discharge.  The  nitrogen  partial  pressure  Is  a 
few  microns  and  the  helium  pressure  is  kept  at  a  few  Torr  to  produce  a 
sufficient  metastable  concentration.  Room  temperature  helium  is  the 
collision  partner  since  the  radiative  decay  curves  from  the  A  ( v ’ =4)  level 
are  sensitive  to  pressure  changes  in  helium  and  not  nitrogen. 

The  probe  laser  scan  of  the  B-X  (5,7)  band  will  resolve  the  nearly 
degenerate  spin  splitting  of  the  X2£g+(vs7,  N)  rotational  levels  be¬ 
cause  the  B 2i;u+(v=5)  manifold  Is  perturbed9*10  by  the  A2^,- (v=17)  level 
of  N2+.  These  perturbations  not  only  cause  shifts  in  the  B  state  rotational 
levels  but  affect  the  relative  intensities  of  the  rotational  structure 
In  the  probe  laser  scan  of  the  B2Eu+-X2£g+(5,7)  band.  The  effects 
of  these  perturbations  on  the  line  strengths  can  be  determined  by 1 1 

J*  1  J  J*  1  J 

S  =  (2J’+1)(2J+1)|Ce  1/2  o  -1/2  +  (uj/ui^e  Cni/2  1/2  -1  1/2  - 


J'  1  J 


Cn3/2  3/2  -1  -1/2  11' 


where  ^  and  ^  are  transition  moments  between  the  X2£g+(v=7)  level  and  the 
A^nu(v=17)  and  B2eu+(v=5)  states,  respectively.  The  3-j  symbols  in 
parentheses  are  evaluated  using  standard  techniques1^.  The  C  coefficients  are 
wavefunction  mixing  factors  that  can  be  detenmined^.lO  from  the  perturbation 
shifts,  spacing  between  the  unperturbed  levels  and  the  non  Bom-Oppenheimer 
matrix  elements  of  the  spin-orbit  {^)  and  orbit-rotation  (n)  interactions, 


K  =  <A2nu|  AL+/2  |B2£g+> 

(2) 

n  =  <A2nu|  Bl+  |B2£g+>  . 

(3) 

The  dominant  factors  in  our  case  is  the  small  value  of  the  ratio  which 

is  approximately1^  o.03  and  the  fact  that  the  rotational  manifolds  do  not 
"cross”9,10  which  result  in  the  intensities14  of  the  perturbed  Rj  branch  being 
reduced  by  only  a  few  percent  for  the  levels  observed  in  our  experiment.  Indeed, 
Klynning  and  Page$9  have  failed  In  their  attempts  to  observe  "extra  lines"  due 
to  this  perturbation.  These  intensity  changes  are  negligible  for  our  present 
experimental  results  but  should  be  considered  in  greater  detail  when  the  signal 
to  noise  ratio  is  substantially  improved  over  our  data. 

Figures  1,  2,  and  3  show  probe  laser  scans  of  the  same  R  branch  portion 
of  the  B-X  (5,  7)  band  with  the  pump  laser  tuned  to  the  Qj(6),  Pi { 8)  and 
02<8)  transitions,  respectively,  of  the  A-X  (4,0)  band.  The  same  helium 
pressure  of  approximately  four  Torr  Is  used  for  these  scans  and  the  pump  and 
probe  laser  pulses  of  about  20  ns  duration  are  coincident  in  time  to  maximize 
single  collision  effects.  The  nearly  degenerate  spin  components  for  the  N 
rotational  levels  of  the  X2Eg+  (v=7)  level  are  clearly  resolved  in  these 


figures  and  the  different  Intensity  patterns  demonstrate  that  propensities 
for  collision  Induced  electronic  energy  transfer  are  strongly  dependent  on 
the  A2^  (v*4,  J)  level  initially  populated.  The  R-|  and  R2  branches 
refer  to  transitions  from  the  (J=N+l/2,  e)  and  F2  (J=N-l/2,  f)  spin 
components,  respectively,  of  the  X2ig+  vibrational  level  and  only  even 
N  values  are  observed  because  of  the  symmetric  (s)  «-*  antisymmetric  (a) 
ruiel.2,8.  The  e  and  f  labels^  refer  to  the  upper  and  lower  series  of 
parity  doublets,  respectively. 
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IV.  Theoretical  Interpretation 

The  propensity  for  aJ  »  0  implies  that  optical  selection  rules  may 
be  important,  but  the  significant  intensities  of  the  higher  |AJ|  lines  in¬ 
dicate  the  presence  of  mechanisms  beyond  dipole  transitions.  Alexander  and  Corey2 
have  developed  a  quantum  mechanical  description  of  inelastic  collision  induced 
transitions  between  2n  and  2E  electronic  states  of  a  diatomic  molecule. 

Their  theory  of  electronic  transfer  by  collisions  is  an  exact  formulation  but 
requires  the  use  of  approximations  to  yield  propensity  rules.  They  use  an 
infinite  order  sudden  ( IOS)  approximation,  which  assumes  a  short  range 
potential  and  collision  energies  large  with  respect  to  the  energy  gaps.  A 
Bom  approximation  may  be  more  appropriate  for  our  experiment  where  the  energy 
gap  is  not  insignificant  compared  to  collision  energies.  The  results  for 
relative  intensity  strengths,  though,  are  similar  in  the  two  limits, 
indicating  that  the  IOS  approximation  may  be  used  to  analyze  our  data. 

Equations  52  and  55a2  for  the  2n  state  in  Hund's  case  (a)  coupling 
may  be  combined  to  yield  the  following  expression  for  the  cross  section  of 
a  transition  between  2n(J)  and  2e+(J')  rotational  states  of  specific 
e ( e/f )  symmetry: 


a  .  *  l  MJ.J'lc 

Je2n  +J'e'2z  A>1  0=t+l/2,G,2n  >  J'=l/2,e.2E 

Q 


where 


k  (J.J*)  =  1/2  [1-  ee'(-l)J+0,+t  ] { 2J ' +1 ) ( 2*+l )(*+!)/( *+2q-1 )  x 


The  3-j  symbols  in  the  bracket  of  Eq.  (5)  determine  the  magnitude  of  the 
coupling  between  the  initial  and  final  rotational  states.  Allowed  values  of 
l  are  determined  by  the  phase  factor  term  in  Eq.  (5)  and  the  triangle  rule  for 
the  magnitudes  of  J,  J'  and  *.  Alexander  and  Corey2  get  propensity  rules  only 
in  the  high-J  limit,  but  predictions  of  relative  line  intensity  strengths  for 
transitions  between  low-lying  rotational  states  may  also  be  obtained  by  evaluating 
the  kA(J,J’)  coefficients  for  specific  rotational  e/f  labels  and  spin  states. 

For  transitions  between  electronic  states  of  opposite  inversion  symmetry,  only  odd 
l  values  contribute  to  the  series  In  Eq.  (1).  Taking  only  the  *=1  term  recovers 
the  dipole  selection  rules,  while  inclusion  of  the  *=3,5,...  terms  adds  necessary 
corrections.  From  Eq.  (4),  the  cross-sections  can  be  determined  only  in  terms  of  the 
*  dependent  "base"  cross-sections,  G( J=*+l/2,e,2n  -*•  J'«1/2,6,2e). 

Thus,  a  priori  cross-sections  cannot  be  obtained  from  this  formulation  but 
relative  values  can  be  determined  only  if  reasonable  estimates  of  the  base  cross 
sections  can  be  extracted  from  the  experimental  data. 

In  Fig.  1,  the  probe  laser's  scan  reflects  the  propensities  for  collision 
induced  transfer  from  the  J=6.5  ( f , s ,-)  level  of  the  A^tu3/2(v=4)  state 
to  the  various  nearly  degenerate  spin  components  of  the  X2£g+(v=7)  ro¬ 
tational  manifold.  Although  the  spin  components  are  clearly  resolved,  the 
separations  are  due  to  perturbations  in  the  B  ^;u+(v=5)  rotational  levels 
so  that  the  relative  intensities  of  a  pair  of  spin  components  in  this  figure 
represent  relative  populations  of  levels  separated  by  ~  0.04  cnr1  in  the 
X2£g+(v=7)  state.  Because  of  the  propensity  for  &J»0  and  our  in¬ 
ability  to  determine  the  contributions  of  the  base  cross  sections  with  *=5 
or  greater,  we  compare  calculated  cross-sections  from  Eqs.  (4)  and  (5)  having 


only  *=1  and  3,  with  the  experimental  data.  In  order  to  get  the  relative 
calculated  cross  sections,  the  ratio  of  the  a(ls D  to  ai Is 3)  base 
cross  sections  must  be  known.  We  chose  this  ratio  to  be  two  In  order  to 
get  the  most  reasonable  agreement  to  the  data.  The  calculated  cross  sections 
predict  that  the  aJ=0,  Rj( 6 .5 )  line  should  be  strongest  with  nearly  twice 
(1.8)  the  intensity  as  the  next  most  Intense  line  which  is  the  aJ=-1 ,  R2< 5.5) 
transition.  Fig.  4(a)  which  compares  the  calculated  and  observed  Intensities 
of  Fig.  1  shows  the  experimental  Rj(6.5)  line  is  slightly  stronger  than 
the  other  lines  but  not  nearly  as  strong  as  predicted.  Most  of  the  other 
observed  lines  are  in  reasonable  agreement  with  the  calculated  Intensities. 

It  should  be  noted  that  In  Figs.  2  and  3,  the  levels  Initially  populated 
In  the  A2nu3/2  and  A2ngl/2  manifolds,  respectively,  have  the  same 
rotational  angular  momentum,  permutation  symmetry  and  parity,  l.e.,  J=7.5(e,s,-) . 
The  Intensity  patterns  for  these  two  figures  are,  however,  quite  different. 

If  we  make  the  same  assumptions  for  the  base  cross-sections  as  In  the  previous 
paragraph,  the  most  Intense  peak  In  Fig.  2  Is  predicted  by  Eqs.  (4)  and  (5)  to  be 
the  aJ=0,  R2< 7.5)  peak  which  Is  expected  to  be  almost  twice  (1.8)  as  large 
as  the  Ri(6.5),  aJ=-1  peak.  In  Fig.  4(b)  the  experimental  intensities  of 
Fig.  2  are  compared  with  their  calculated  values.  As  In  Fig.  4(a),  the 
experimental  R2 (7.5)  peak  In  Fig.  4(b)  does  not  agree  with  its  predicted 
Intensity.  The  other  peaks,  however,  are  within  15%  of  the  predicted  relative 
cross  sections  which  is  excellent  agreement.  Fig.  4(c)  gives  the  comparison 
of  calculated  and  observed  Intensities  of  Fig.  3.  The  strongest  predicted 
line  Is  the  a«J=0,  R2( 7.5)  transition  with  nearly  twice  the  Intensity  as 
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the  aJ=+1,  8.5)  line  which  should  have  the  next  largest  predicted 

cross  section.  The  figure  shows  that  R2 ( 7 - 5 >  is  Indeed  the  strongest 
line,  but  the  next  most  Intense  line  Is  not  Ri(8.5).  For  this  particular 
case  where  a  strong  line  occurs  as  predicted,  the  weaker  lines  are  not 
In  as  good  agreement  as  the  above  cases. 

In  general,  probe  laser  scans  of  collision  Induced  transitions  from 
rotational  levels  having  the  same  e  or  f  label  in  a  particular  spin  component 
of  the  state  have  similar  intensity  patterns.  Thus,  the  comparisons 

given  above  are  representative  of  the  degree  of  agreement  between  our  experi¬ 
mental  results  and  the  relative  cross  sections  based  on  the  work  of  Alexander 
and  Corey2.  Considering  the  assumptions  that  had  to  be  made  concerning  the 
base  cross  sections  In  order  to  calculate  the  relative  cross  sections  from 
Eqs.  (4)  and  (5),  we  find  our  results  to  be  In  general  qualitative  agreement 
with  the  predicted  Intensities.  There  Is  still  need,  however,  for  more 
theoretical  as  well  as  experimental  work  before  these  collision  Induced 
electronic  transition  processes  can  be  understood.  It  Is  still  not 
clear  how  these  electronic  transfers  can  occur  over  a  large  energy  gap 
(  -1760  cm-1)  with  the  efficiency  necessary  to  reveal  propensity  rules. 

At  present,  work  on  determlng  the  form  of  the  Interaction  potential 
Is  continuing. 
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V.  Recommendations 


Obviously,  the  work  started  in  this  research  has  yet  to  be  completed. 
As  has  been  seen  in  the  previous  sections,  the  differences  between  theory 
and  experiment  are  still  great.  A  simple  quantum  mechanical  formulation  such 
as  that  of  Alexander  and  Corey,  though  intuitive  and  complete,  still  lacks 
the  ability  to  quantitatively  predict  the  relative  cross  sections  for  the 
transitions  observed  in  our  experiment. 

The  inclusion  of  the  interaction  potential  into  these  calculations  seems 
to  be  essential.  Further  work  by  others  has  shown  that  pursuing  the  close 
coupling  formulation  without  specifying  the  potential  leads  to  no  better 
results  than  that  oDtained  in  considering  only  the  first  order  approximation 
of  the  theory. 

Thus  far,  we  have  not  been  able  to  arrive  at  theoretical  results  which 
will  support  our  experimental  findings,  but  research  continues  in  this  area. 
Other  improvements  in  the  theory  are  possible  that  may  alter  the  existing 
predictions.  One  example  would  be  to  examine  a  more  involved  trajectory 
of  the  scattering  partner  in  the  course  of  one  of  these  collisions.  Up  until 
the  present,  a  straight-line  approximation  has  been  used  to  make  the  cal¬ 
culations  time-dependent  and  thereby  include  the  effects  of  the  energy  gap; 
however  such  an  assumption  is  by  no  means  realistic. 
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Figure  Captions 


Fig.  1  OODR  spectrum  of  the  B +-X  ^g  +  (5,7)  band  of  N2+. 

The  numbers  of  the  Rj  and  R2  branches  refer  to  the  J"  levels 
of  this  band.  The  probe  laser  scans  this  band  with  the  pump 
laser  tuned  to  the  Qt(N,,=6,  J"=6.5)  line  of  the  A-X  (4,0) 
band  so  that  the  A^nu3/2(v=4,  j=6.5,  f,  s,  -)  level 
is  selectively  populated. 

Fig.  2  OODR  spectrum  of  the  B-X  (5,7)  band  with  the  pump  laser  tuned 
to  the  Pi(N"=8,  J"=8.5)  line  of  the  A-X  (4,0)  band.  The 
A2nu3/2^v=4*  0*7.5,  e,  s,  -)  level  is  initially  populated 
for  this  scan. 

Fig.  3.  OODR  spectrum  of  the  B-X  (5,7)  band  with  the  A^u-|y2^v=A» 

J=7.5,  e,  s,  -)  level  being  selectively  populated  by  the  pump 
laser  tuned  to  the  O2  (N"=8,  J"=7.5)  line  of  the  A-X  (4,0)  band. 

Fig.  4.  Comparison  of  theoretical  to  experimental  intensities  obtained  with 
the  pump  laser  tuned  to  (a)  Qi ( 6 ) ,  (b)  P ^ ( 8 ) ,  and  (c)  02(8)  lines 

O  O  1 

of  the  A^riy-X^Eg  (4,0)  band.  The  corresponding  probe  laser 
scans  which  show  these  experimental  intensities  are  in  Figs.  1, 

2,  and  3,  respectively.  For  each  R( N" )  spin  doublet,  Rj( J"=N"+l/2) 
lies  to  the  left  of  the  R2( J"=N"-l/2)  component. 
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In  the  one  month  period  -following  the  completion  of  my 
first  vear  of  research,  the  advances  have  been  completely  in  the 
theoretical  modelling  of  our  collision-induced  electronic  transition 
experiment.  As  was  previously  reported.  Dr.  Daniel  H.  Katayama 
and  I  have  completed  an  opt i cal -opt i c al  double  resonance  experiment 
in  which  we  were  able  to  observe  transitions  between  specific 
rovi  br  at  i  onal  manifolds  in  the  electronic  transition  AiTTw-Xx£^’  in 
collisions  of  N*  with  atomic  helium.  No  previously  suggested  model 
has  been  able  to  reproduce  our  data,  and  so  we  have  begun  trying 
to  model  the  interaction  ourselves. 

Assuming  a  Morse-type  interaction  potential  similar  to  that 
currently  used  in  describing  vibrational  pred i ssoc i at i on  of  diatomic 
molecules,  we  have  developed  an  electronic  analogue  to  compare 
to  our  experimental  findings.  With  this  particular  model,  we  will 
be  able  to  study  the  sensitivity  of  our  proposed  interaction  to 
collision  angle,  relative  orbital  angular  momentum  before  and 
after  collision  and  energy  gap  traversed  in  transition. 
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Neutral  Reactions  in  the  Presence  of  Alkali  Ions 


A.A.  Viggiano1  ,  Carol  A.  Deakyne2,  F.  Dale,  and  John  F.  Paulson 
AFGULID  Hanscom  AFB  MA  01 731  -5000 

ABSTRACT 

It  has  been  shown  that  neutral  ligands  clustered  to  alkali  ions  can  react  with 
other  neutrals.  In  fact,  the  presence  of  the  alkali  ion  greatly  enhances  the  rate  of  the 
corresponding  neutral  -  neutral  reaction.  Several  more  reactions  that  exhibit  this 
phenomenon  have  been  identified  experimentally,  with  rate  enhancements  ranging 
from  4  to  30  orders  of  magnitude.  Three  possible  explanations  for  the  rate 
enhancement  have  been  explored.  Ab  initio  calculations  have  been  carried  out  to 
supplement  the  experimental  measurements  and  to  provide  further  insight  into  the 
mechanism.  The  calculations  indicate  that  bonding  the  neutral  to  the  alkali  ion  causes 
structural  rearrangements  that  shift  the  geometry  of  the  neutral  molecule  toward  its 
transition  state  geometry,  thereby  lowering  the  activation  energy  of  the  reaction. 


1.  Under  contract  to  AFGLfrom  Systems  Integration  Engineering  Inc.,  Lexington  MA 
02173. 

2.  Air  Force  Geophysics  Scholar 


Introduction 

Alkali  ions  are  isoelectronic  with  the  noble  gases  and  therefore  are  unreactive 
species.  However,  they  form  cluster  bonds  to  a  wide  variety  of  neutrals.  These  cluster 
bonds  have  been  shown  to  be  mainly  electrostatic  in  nature.1  Rowe,  et  al.2  found  that 
the  ligands  clustered  to  the  alkali  ions  can  react  with  other  neutrals,  and,  in  fact,  the 
presence  of  the  alkali  ion  greatly  enhances  the  rate  of  the  corresponding  neutral  - 
neutral  reaction.  For  example,  the  rate  constant  for  reaction  (1)  is  9  orders  of 
magnitude  larger  in  the  presence  of  a  Li"*'  ion  than  the  corresponding  neutral  reaction.2 

(1)  N205  +  N0  3N02 

The  rate  enhancement  for  this  reaction  was  found  to  depend  strongly  on  the  alkali  to 
which  the  N205  was  clustered.  The  rate  constant  for  reaction  (2)  was  also  significantly 
enhanced  in  the  presence  of  alkali  ions. 

(2)  Oj  +  NO  — >  N02  +  02 

Rowe,  et  al.2  postulated  three  reasons  for  these  large  rate  enhancements.  The 
first  is  that  the  electrostatic  interaction  energy  of  the  alkali  ion  -  neutral  complex  and 
reactant  neutral  is  larger  than  the  activation  energy  for  the  reaction  between  the  neutral 
reactants  in  the  absence  of  the  ion.  The  second  is  that  the  potential  energy  surface  is 
sufficiently  altered  in  the  presence  of  the  ion  to  greatly  increase  the  rate  constant. 
Finally,  the  lifetime  of  the  collision  complex  in  the  presence  of  the  alkali  ion  is  longer, 
yielding  a  greater  reaction  probability. 

We  have  untertaken  a  study  to  determine  the  relative  importance  of  the  above 
causes  of  the  rate  enhancement.  We  have  found  several  more  reactions  that  exhibit 
this  phenomena,  with  rate  enhancements  of  as  much  as  30  orders  of  magnitude.  In 
addition,  ab  initio  calculations  have  been  carried  out  in  order  to  supplement  the  rate 
constant  measurements  and  to  provide  further  insight  into  the  mechanism. 


A.  Experimental  Method.  The  experimental  measurements  were  made  in  a 
newly  constructed  combination  flowing  afterglow  (FA)  -  selected  ion  flow  tube  (SIFT) 
apparatus.  The  apparatus  is  shown  in  Figure  1.  The  SIFT  part  of  the  instument  is 
similar  to  the  previous  instrument  used  in  our  laboratory.  The  FA  was  designed  to 
operate  at  high  pressures  using  a  design  similar  to  that  of  Fahey  et  al.3  In  the  SIFT, 
ions  are  created  in  a  high  pressure  (-1  torr)  ion  source.  Upon  exiting  the  ion  source, 
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the  ions  are  focused  into  a  quadrupole  mass  filter  and  are  injected  through  a  venturi 
inlet  into  a  flow  tube  1  meter  long.  Introducing  the  ions  through  the  venturi  inlet  aids 
the  movement  of  ions  from  the  low  pressure  quadrupole  region  to  the  high  pressure 
flow  tube.  The  ions  are  carried  down  the  length  of  the  flow  tube  by  the  buffer  gas.  The 
neutral  reactant  is  added  downstream.  Reactant  and  product  ions  are  sampled  though 
a  0.2  mm  hole  in  a  blunt  nose  cone  and  detected  by  a  second  quadrupole  and  channel 
electron  multiplier. 

The  FA  differs  only  in  that  the  ions  are  created  directly  in  the  carrier  gas.  For 
these  studies  the  ion  source  region  was  separated  from  the  flow  tube  by  a  diaphragm  5 
mm  in  diameter.  Under  these  conditions  the  ion  source  region  was  maintained  at  a 
pressure  of  10-20  torr  and  the  flow  tube  was  operated  at  approximately  0.5  torr. 

The  alkali  ions  were  produced  by  thermionic  emission  from  a  rhenium  filament 
coated  with  the  appropriate  alkali  nitrate,  silicon  dioxide  and  aluminum  oxide.  After  a 
few  hours  of  conditioning,  heating  the  filament  produced  an  essentially  pure  signal  of 
the  alkali  of  interest.  The  filament  was  electrically  biased  by  a  few  volts  with  respect  to 
the  flow  tube  walls  in  order  to  help  the  ions  enter  the  flow.  Three  filaments  were  put  in 
the  tube  to  enable  us  to  change  alkali  ions  rapidly.  The  cluster  ions  were  made  by 
adding  the  clustering  neutral  to  the  ion  source  region  slightly  downstream  of  the 
filament.  For  ihe  most  part,  clustering  occurred  in  the  ion  source  region  due  to  the  high 
pressure  there.  In  this  way,  we  were  frequently  able  to  add  much  less  of  the  clustering 
gas  than  would  have  been  needed  to  obtain  the  same  amount  of  clustering  in  the 
absence  of  the  diaphragm.  In  order  to  obtain  the  best  signal  for  the  cluster  ion  of 
interest  without  interference  from  other  ions  (i.e.  larger  cluster  ions),  the  temperature, 
pressure,  and  buffer  gas  were  varied.  Three  buffers  were  used,  He,  N2,  and  Ar. 

The  gases  used  in  the  experiments  were  standard  commercial  gases.  The  CO 
and  NO  were  purified  by  passing  them  through  an  ascarite  filter.  Atomic  oxygen  was 
produced  by  passing  a  mixture  of  He  and  02  through  a  microwave  discharge.  The 
amount  of  O  formed  was  monitored  by  the  reaction  of  CH5+  (generated  in  the  SIFT) 
with  O.4  Typically  the  percent  dissociation  was  in  the  30-40%  range.  This  method  is 
not  as  accurate  as  measuring  the  O  atom  concentration  directly  but  is  sufficient  for  the 
present  purposes.  Rate  constants  are  believed  to  be  accurate  to  within  a  factor  of  two. 

B.  Theoretical  Method.  The  calculations  were  carried  out  ab  initio  on  a  VAX 
11/780  computer.  The  Gaussian  82  series  of  programs  were  utilized  to  perform  the 
calculations.5  The  optimum  structures  of  C02,  N02,  NO,  Li+(C02).  Li+(N02),  Li+(NO), 
and  LiO+  were  obtained  at  the  HF/6-31G*  basis  set  level6  via  the  force  relaxation 
method.7  All  the  additional  calculations  were  carried  out  at  the  optimized  6-31 G* 
geometries.  Electronic  lithium  ion  affinities  AE's  were  computed  with  the  6-31 G*  and 
6-31 +G*  basis  sets.8  The  effect  of  electron  correlation  on  the  AE  values  was 
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determined  by  means  of  Metier  -  Plesset  perturbation  theory  to  third  (MP3)  order.9 

From  her  investigation  of  the  lithium  ion  affinities  of  several  oxygen  and  nitrogen 
bases,  Del  Bene10  has  shown  that  relative  affinities  at  these  levels  of  calculation  are 
quite  good.  However,  the  AE  values  tend  to  be  overestimated  (particularly  the  6-31 G* 
AE's)  compared  to  the  values  obtained  at  higher  levels  of  theory. 

The  electronic  lithium  ion  affinities  were  calculated  using  equation  (3).  Here  ET 
is  the  total  electronic  energy  and  A  is  0,  NO,  CO,  N02,  or  C02.  In  order  to  compare  the 
AE’s  to  experimental  gas  -  phase  lithium  ion  affinities,  zero  -  point  vibrational  energy 
changes  and  temperature  effects11  have  been  estimated  based  on  Del  Bene’s10 
results. 


(3)  AE  » -  [Et  (LiA+)  -  Et  (Li+)  -  ET  (A)] 

Results  and  Discussion 

A.  Trends  in  the  Experimental  Rate  Enhancements.  The  results  of  this  study  as 
well  as  those  of  Rowe  et  al.2are  listed  in  Table  I.  In  addition  to  the  reactions  given  in 
the  table  Rowe  et  al.2  studied  two  other  reactions  that  showed  no  rate  enhancement. 
None  of  the  rate  constants  were  found  to  vary  significantly  with  temperature.  The 
column  in  Table  I  headed  rate  enhancement  refers  to  the  ratio  of  the  rate  constant  in 
the  presence  of  the  alkali  ion  to  that  of  the  corresponding  neutral  reaction  at  the  coldest 
temperature  studied.  Since  none  of  the  rate  constants  reported  in  Table  I  have  been 
found  to  vary  with  temperature  and  all  of  the  neutral  rate  constants  have  positive 
activation  energies,  this  definiton  represents  the  largest  value  of  the  rate  enhancement 
for  each  reaction.  The  rate  enhancement  varies  from  4  orders  of  magnitude  to  30 
orders  of  magnitude. 

As  a  result  of  several  different  experimental  difficulties  we  were  unable  to 
measure  any  of  these  reactions  on  K+.  In  general,  bond  strengths  of  neutrals  to  alkali 
ions  decrease  as  the  molecular  weight  of  the  alkali  ion  increases,  due  to  the  increasing 
size  of  the  alkali.12  Therefore,  one  must  decrease  the  temperature  of  the  flow  tube  in 
order  to  get  workable  signals  of  the  K+  clusters.  In  the  case  of  K+(N02)  the  N02  froze 
in  the  inlet  line  before  we  could  obtain  usable  signals  of  K+(N02). 

We  could  not  do  any  experiments  that  involved  K+  and  O  due  to  the  following 
unusual  problem.  Some  unknown  neutral  KX  emitted  from  the  K  filament  reacted  with 
O  to  form  ion  pairs. 


We  could  detect  the  formation  of  both  K+  and  X(0)\  Negative  ion  masses  were 
observed  at  79±0  and  246 ±2  amu  in  the  ratio  of  1  to  10.  We  have  not  been  able  to 
identify  these  ions  although  the  ion  at  mass  79  has  no  major  isotopes.  Reaction  (4) 
has  a  rate  constant  of  2x1  O'10  cm3  s'1  and  unfortunately  produces  much  more  K+  than 
is  emitted  directly  from  the  filament.  We  had  also  hoped  to  study  several  reactions 
involving  F2  but  neutrals  emitted  from  all  three  alkali  filaments  created  ion  pairs  upon 
addition  of  F2.  These  chemiionization  reactions,  although  not  understood,  represent 
the  majority  of  the  ground  state  chemiionization  reactions  known.13  Moreover,  they 
are  the  only  chemiionization  reactions  that  form  negative  ions  directly  rather  than 
forming  electrons.13 

Several  trends  in  the  data  are  readily  observable.  Neutrals  clustered  to  Li+ 
react  either  at  the  same  rate  or  faster  than  those  clustered  to  Na+,  which  in  turn  react 
faster  than  those  clustered  to  K+  (Table  I).  The  most  dramatic  example  is  the  case  of 
N205  reacting  with  NO,  where  the  reaction  in  the  presence  of  Li+  is  about  two  orders  of 
magnitude  faster  than  it  is  in  the  presence  of  Na+.  No  detectable  reaction  was 
observed  in  the  presence  of  K+.  The  rate  constant  observed  for  the  the  reaction  of  N02 
with  CO  in  the  presence  of  Li+  is  a  factor  of  seven  slower  than  it  is  in  the  presence  of 
Na+.  For  the  03  reaction  with  NO  the  rates  in  the  presence  of  Li+  and  Na+  are  the 
same,  while  the  K+  reaction  is  over  an  order  of  magnitude  slower.  All  the  reactions 
involving  O  atoms  are  fast  (within  a  factor  of  5  of  the  collison  rate)  with  the  exception  of 
the  reaction  of  N20  with  O. 

In  some  cases  the  two  neutral  products  both  separate  from  the  alkali  as,  for 
example,  in  the  03  reaction  with  NO.  This  reaction  is  highly  exothermic  which  provides 
sufficient  energy  to  dissociate  both  neutrals  and,  indeed,  both  do  dissociate.  The 
reaction  of  N02with  CO  is  also  sufficiently  exothermic  to  dissociate  both  ligands; 
however,  in  this  case  the  C02  ligand  stays  attached  to  the  alkali  ion.  In  the  reaction  of 
N205  with  NO  the  reaction  is  exothermic  only  if  one  of  the  products  remains  attached  to 
the  alkali  ion,  and  this  is  what  is  observed. 

For  several  of  the  reactions  it  was  also  possible  to  determine  what  occurs  when 
a  second  ligand  is  bonded  to  the  alkali  ion  -  neutral  complex  (Table  I).  The  rate 
constant  was  found  to  be  approximately  the  same  whether  one  or  two  ligands  were 
bonded  to  the  alkali  ion,  although  the  product  of  the  two  reactions  sometimes  varied. 
For  example,  when  H2S  reacts  with  O  the  main  products  are  Li+(HSO)  (or  Na+(HSO)] 
+  H  when  one  ligand  is  attached  and  Li+(H2SO)  [or  Na+(H2SO)J  +  H2S  when  two 
ligands  are  attached.  The  reaction  involving  two  ligands  may  be  thought  of  as  an 
addition  reaction  in  which  the  third  body  (the  second  H2S)  is  built  into  the  reactants. 

B.  Trends  in  the  Calculated  Electronic  Lithium  Ion  Affinities.  Table  II  tabulates 
the  total  energies  of  Li+,  O,  NO,  N02,  C02,  LiO+,  Li+(NO),  Li+(N02),  and  Li+(C02)  at 
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the  various  basis  set  levels.  Table  III  gives  the  lithium  ion  affinities  of  0,  NO,  CO,10 
N02,  and  C02.  LiO+  is  essentially  not  bound  at  any  level  of  theory.  The  relative 
affinities  of  the  other  species  are  C02  >  N02  >  CO  >  NO. 

The  addition  of  diffuse  functions  to  the  basis  set  (6-31 +G*  vs.  6-31 G,  Table  III) 
decreases  AE  by  about  1-1.5  kcal  mol'1.  Taking  electron  correlation  effects  into 
account  at  the  MP2  level  also  decreases  AE  by  about  2  -  3.5  kcal  mol'1.  In  contrast,  the 
third  order  Maller-Plesset  calculation  increases  the  lithium  ion  affinity  to  a  value  in 
between  the  Hartree  -  Fock  (HF)  and  MP2  values.  These  patterns  are  identical  to 
those  observed  by  Del  Bene.10 

The  MP3/6-31+G*  electronic  affinities  obtained  by  Del  Bene10  are  1-1.5  kcal 
mol'1  larger  than  those  obtained  with  the  highest  level  of  theory  she  considered,  i.e. , 
the  MP4SDQ/6-31 1  +G(2d,2p)  level.  In  addition,  correcting  for  zero  -  point  energy 
changes  and  temperature  effects  (i.e.,  converting  a  AE  to  a  AH)11  lowers  Del  Bene's 
calculated  AE's  by  0.5  -  1.5  kcal  mol*1.  If  the  lithium  ion  affinities  of  the  molecules 
studied  in  this  work  follow  the  same  trends  as  those  observed  by  Del  Bene,10  then 
increasing  the  size  of  the  basis  set,  taking  the  Mailer  -  Plesset  calculations  to  higher 
order,  and  correcting  for  zero  -  point  energy  changes  and  temperature  effects  will 
decrease  the  values  reported  in  Table  II  by  1  -  3  kcal  mol'1.  Corrected  lithium  ion 
affinities  computed  with  the  MP4SDQ/6-31  t+G(2d,2p)  basis  set  are  generally  about  2 
kcal  mol'1  smaller  than  the  experimental  gas  -  phase  values.14 

Although  changing  the  basis  set  does  change  the  magnitudes  of  the  AE's,  the 
trends  in  the  AE's  vary  very  little  with  respect  to  basis  set  (Table  III).  The  order  of  the 
electronic  lithium  ion  affinities  is  the  same  at  each  basis  set  level,  and  the  5AE’s  are 
quite  similar  at  each  basis  set  level. 

C.  Reaction  Products  vs.  Buffer.  One  interesting  observation  not  included  in 
Table  I  is  that  the  products  of  the  N02  reaction  with  CO  in  the  presence  of  Li+  vary 
depending  on  which  buffer  is  used.  In  a  He  buffer  the  major  product  ion  is  Li+(CO);  the 
minor  product  ion  is  Li+(C02).  In  a  N2  buffer  only  the  Li+(C02)  product  is  obtained. 
About  equal  amounts  of  the  two  products  are  formed  in  an  Ar  buffer. 

The  following  explanation  is  consistent  with  the  above  observations,  although 
we  have  no  firm  evidence  that  excludes  other  explanations.  Reaction  (5) 

(5)  U+(N02)  +  CO  -»  Li+(C02)  +  NO 

is  highly  exothermic;  AH  is  estimated  to  be  -73  kcal  mol'1  from  the  HF/6-31G*  data  in 
Table  II  [ET(CO)  =  -112.73691  kcal  mol'1  15J.  This  suggests  that  the  initial  products 
contain  a  considerable  amount  of  internal  energy.  Helium  is  not  a  very  good  quencher 
of  vibrational  energy  in  ions;  therefore,  one  would  expect  the  Li+(C02)  product  to 
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remain  vibrationally  excited.16  In  contrast,  one  would  expect  that  any  vibrational 
excitation  would  be  quenched  in  a  N2  buffer  and  that  an  Ar  buffer  would  quench  some 
but  not  all  of  the  excitation.16  The  switching  reaction 

(6)  U+(C02)  +  CO  ->  Li+(CO)  +  C02 

is  endothermic.  AH  is  estimated  to  be  3  -  5  kcal  mol'1  based  on  the  data  in  Table  III 
and  the  analysis  given  above  for  converting  AE  to  AH.  However,  if  the  reactant  ion  has 
vibrational  energy  in  excess  of  the  endothermicity,  the  reaction  may  proceed.  In  light  of 
the  above  discussion  on  quenching  rates  in  different  buffers,  one  would  expect  the 
amount  of  Li+(CO)  formed  to  vary  in  the  manner  observed. 

D.  Reasons  for  .the  Rate  Enhancement.  One  would  like  to  obtain  some 
understanding  of  the  reasons  for  the  rate  enhancement.  Only  a  detailed  theoretical 
study  can  answer  this  interesting  problem  but  some  simple  qualitative  arguments  may 
shed  some  light  on  it.  Three  explanations  have  been  postulated  for  why  the  reaction 
rate  constants  are  significantly  faster  in  the  presence  of  the  alkali  ions.2  The  first  is  that 
the  interaction  energy  between  the  cluster  ion  and  neutral  is  significantly  larger  than 
that  for  two  neutrals.  Secondly,  the  moiecular  potential  energy  surface  is  signifantly 
altered  by  the  presence  of  the  ion.  Finally,  the  lifetime  of  an  ion  -  neutral  collision 
complex  is  significantly  longer  than  that  for  a  neutral  -  neutral  encounter.  In  reality 
these  are  not  separable  quantities  but  it  is  of  interest  to  look  at  the  information  we  have 
on  each  one  individually. 

The  interaction  energies  between  the  ionic  and  neutral  reactants  are  not  known. 
Two  simple  methods  are  available  to  estimate  the  order  of  these  energies.  The  first 
method  is  to  take  the  interaction  energy  between  a  positive  charge  and  a  polarizable 
dipole  at  the  radius  of  the  alkali  ion.  Assuming  purely  electrostatic  interactions  this 
would  overestimate  the  interaction  energy,  since  the  clustered  alkali  has  a  larger 
distance  of  closest  approach  than  the  unclustered  alkali.  The  other  method  is  to  use 
information  on  the  bond  strengths  of  other  ligand  -  alkali  ion  complexes  to  obtain  some 
estimate  of  the  bond  strengths  of  the  species  in  question.  The  first  method17  gives 
interaction  energies  of  45  -  129  kcal  mol'1  for  Li+,  7  -  26  kcal  mol'1  for  Na+  and  2  -  9 
kcal  mol'1  for  K+  with  the  neutrals  used  in  this  study.  The  lower  number  refers  to  0 
atoms  and  the  higher  to  CO.  Not  much  information  is  available  on  bond  strengths  of 
neutrals  to  alkali  ions.  Typical  bond  strengths  of  Li+  to  polar  neutrals  are  on  the  order 
of  35  kcal  mol'1.12  Na+  clusters  have  bond  strengths  between  12  (CO)  and  24  (H20) 
kcal  mol'1.)12  K+  bond  strengths  range  from  8.5  (C02)  to  17  (H20)  kcal  mol'1.12 

A  first  estimate  of  whether  the  neutral  reaction  can  take  place  in  the  presence  of 
the  alkali  ions  is  whether  the  above  electrostatic  interaction  energies  are  larger  than 
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the  neutral  activation  energies  Ea.  In  fact,  this  was  the  method  used  to  determine 
which  reactions  might  be  catalyzed  by  the  presence  of  the  alkali  ions.  Previously,  no 
such  method  was  employed  and  most  of  the  reactions  tried  were  not  catalysed. 

The  activation  energy18  for  the  reaction  between  N205  and  NO  is  22  kcal  mol'1. 
The  estimated  Alk+(N205)  -  NO  interaction  energies  are  >  35  kcal  mol'1  for  Li+,  15-25 
kcal  mol'1  for  Na+  and  £10  kcal  mol'1  for  K+.  These  estimates  are  in  line  with  the 
observed  trends.  Li+  enhances  the  rate  the  most,  and  the  Li+  interaction  energy  is 
much  larger  than  the  activation  energy.  The  rate  constant  in  the  presence  of  Na+  is 
two  orders  of  magnitude  lower  than  it  is  in  the  presence  of  Li+,  and  the  interaction 
energy  and  activation  energy  are  comparable  for  Na+.  The  rate  is  not  affected  in  any 
measurable  manner  when  K+(N205)  rather  than  N205  reacts  with  NO.  In  this  case  the 
interaction  energy  is  significantly  lower  than  the  activation  energy. 

The  rate  constants  for  the  reactions  of  Na+(03)  and  Li+(03)  with  NO  are  similar 
in  magnitude  (Table  I),  while  that  for  the  reaction  of  K+(03)  with  NO  is  an  order  of 
magnitude  smaller.  The  activation  energy19  for  the  neutral  reaction  between  03  and 
NO  is  2.9  kcal  mol'1.  The  electrostic  interaction  energies  are  significantly  larger  than 
Ea  for  Li+  and  Na+.  In  contrast,  Ea  and  the  K+(03)  -  NO  interaction  energy  are  similar 
(see  analysis  above).  For  the  0  reactions  (except  for  N20  +  O),  all  of  which  are  fast, 
the  interaction  energies  are  significantly  greater  than  the  activation  energies  of  3  -  4 
kcal  mol'1.20  The  reaction  of  N02  with  CO  has  the  largest  activation  energy  of  any 
reaction  for  which  we  observed  a  rate  enhancement.  Here  Ea  is  less  than  the 
Alk+(N02)  interaction  energy  expected  for  Li+  but  is  slightly  larger  than  that  expected 
for  Na+.  We  also  looked  at  the  reaction  of  H20  +  CO  which  has  an  activation  energy  of 
53  kcal  mol'1  and  saw  no  reaction  in  the  presence  of  Ll+.  Rowe  et  al.2  also  studied 
several  reactions  with  large  activation  energies  and  saw  no  reaction. 

E.  Reactant  Surfaces.  The  above  results  may  best  be  explained  by  a  change  in 
the  reactant  surface  when  the  neutrals  are  bound  to  an  alkali  ion  (Figure  2).  The  top 
drawing  in  Figure  2  is  a  schematic  representation  of  the  reaction  coordinate  for  the 
purely  neutral  reaction.  Reactants  A  and  B  proceed  over  an  activation  barrier  of 
magnitude  Ea  yielding  reactants  C  and  D.  The  simplest  model  for  the  reaction  in  the 
presence  of  the  alkali  ion  is  the  two  basin  potential  proposed  by  Brauman  and 
co-workers21  to  explain  many  ion-molecule  reactions  (bottom  drawing,  Figure  2).  The 
reactants  are  neutral  A  clustered  to  an  alkali  ion  Alk+(A)  and  B.  As  the  reaction 
proceeds  along  the  reaction  coordinate  the  reactants  fall  into  an  attractive  well  due  to 
the  electrostatic  attraction  between  the  cation  cluster  and  neutral.  After  the  first 
attractive  well,  Alk+(A)(B)  encounters  a  barrier  whose  height  can  be  either  above  or 
below  the  initial  energy  of  the  reactants.  A  second  well  follows  corresponding  to  the 
cluster  Alk+(C)(D).  Finally  one  or  both  of  the  neutral  products  dissociates  form  the 
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alkali  ion. 

Based  on  the  above  schematic  the  criterion  for  reaction  is  the  height  of  the 
activation  barrier  between  the  two  wells.  If  the  barrier  exceeds  the  initial  energy  of  the 
reactants  (case  A)  then  the  reaction  will  still  have  a  positive  activation  energy  and  will 
not  be  measurable  in  our  apparatus.  The  reaction  may  still  be  significantly  faster  in  the 
presence  of  the  ion  but  we  cannot  measure  the  difference  in  rates.  If  the  barrier  is 
considerably  lower  than  the  initial  energy  of  the  reactants  (case  B)  the  reaction  rate 
should  be  rapid.  As  the  barrier  height  changes  from  above  to  below  the  initial  energy 
of  the  reactants  (case  C)  the  rate  constant  should  become  increasingly  larger.  In 
nucleophilic  displacement  reactions  involving  gas  phase  ions  the  effect  of  the  barrier 
becomes  negligible  only  when  it  is  approximately  5-10  kcal  mol'1  below  the  zero  of 
energy.21 

How  does  the  present  data  fit  into  this  scheme?  All  the  reactions  in  the 
presence  of  Li+  would  be  represented  by  case  B,  i.e. ,  the  barrier  has  little  or  no  effect 
since  the  well  depth  is  substantial.  The  reactions  in  the  presence  of  Na+  for  which  the 
rate  constant  is  comparable  to  that  for  the  reaction  in  the  presence  of  Li+  are  also  type 
B  reactions.  Those  for  which  the  rate  is  substantially  slower  than  the  rate  when  a  Li+ 
complex  is  a  reactant  would  be  case  C  reactions.  The  barrier  height  has  just  begun  to 
be  important  when  the  two  rate  constants  vary  by  less  than  an  order  of  magnitude.  The 
reactions  of  N205  and  03  with  NO  in  the  presence  of  K+  fall  into  categories  A  and  B, 
respectively.  As  stated  previously  experimental  difficulties  prevented  us  from  studying 
any  of  the  new  reactions  in  the  presence  of  K+.  This  is  unfortunate  since  it  would  have 
been  a  useful  test  of  this  scheme. 

The  temperature  dependences  of  reactions  in  the  three  categories  would  be  as 
follows.  Case  A  reactions  would  have  a  positive  temperature  dependence,  case  B 
reactions  would  show  no  temperature  dependence,  and  case  C  reactions  would  have 
a  negative  temperature  dependence.  The  magnitude  of  the  the  negative  temperature 
dependence  would  be  fairly  small  ranging  between  approximately  T",/2  (small 
molecules  only)  and  T°,  where  I  is  the  total  number  of  rotational  degrees  of  freedom.22 
No  temperature  dependences  were  seen  for  any  of  the  reactions,  although  a  small 
temperature  dependence  would  have  been  missed  due  to  the  larger  than  normal 
uncertainty  in  the  data.  This  is  in  qualitative  agreement  with  the  above  picture. 

L _ vs.  E^'.  Next  we  explore  the  interesting  problem  of  the  relative 

magnitudes  of  the  barrier  heights  Ea  and  Ea'.  The  discussion  above  correlating  the 
observed  rate  constant  to  the  electrostatic  interaction  energy  and  the  neutral  activation 
energy  suggests  that  Ea  and  Ea'  are  similar  but  not  identical  in  magnitude. 

Ab  Initio  Structures.  In  order  to  gain  some  insight  into  how  the  activation  barrier 
is  affected  by  the  presence  of  the  ion  consider  the  ab  initio  calculations  on  the  structure 
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and  energetics  of  the  complexes.  The  Li+(NO)  structure  was  optimized  in  the  Cs  point 
group  and  remained  bent  with  ZLi  •••  O  -  N  =  175.2°.  The  Li  •••  O  distance  is  1 .965  A 
and  the  O  -  N  distance  is  1.137  A.  The  latter  distance  can  be  compared  with  the  bond 
length  of  1.127  A  in  neutral  NO.  Li+(C02)  was  optimized  in  the  Cs  and  C2v  point 
groups.  The  bent  Cs  structure  is  not  a  local  minimum;  it  changes  to  the  linear  C2v 
structure  upon  optimization. 

Four  possible  conformations  of  Li+(N02)  were  examined  and  are  shown  in  the 
diagram  below.  Their  point  groups  are  also  given  in  the  diagram.  The  most  stable 
structure  is  II.  The  equilibrium  geometries  of  C02,23  N02,  Li+(C02),  and  form  II  of 
Li+(N02)  are  presented  in  Figure  3.  Notice  that  there  is  a  correlation  between  the 
electronic  lithium  ion  affinities  and  the  Li  •••  0  bond  lengths.  As  the  Li  •••  O  bond  length 
decreases,  AE  increases  (Table  III  and  Figure  3). 


U 

o 

X 

nJ 

0  y 

Ll  -  °"  KJ 

Cav 

n,  cs 

Ll  •’ 

■ 

ro 

i 

Ll  *•' 

£LT 

,Cs 

C 

JV}  c, 

Compare  the  structures  of  N02  and  Li+(N02)  in  Figure  3.  In  order  to  facilitate 
the  reaction  N02  +  CO  -»  C02  +  NO,  it  is  desirable  to  lengthen  (i.e.,  weaken)  one  of  the 
NO  bonds  and  to  shorten  the  other.  That  is  precisely  what  occurs  in  the  presence  of 
the  Li+  ion.  In  fact,  the  same  type  of  structural  rearrangement  occurs  in  the  Li+(C02) 
cluster  (Figure  3)  and  in  other  lithium  clusters,  also.24*25  This  suggests  that  forming 
the  Alk+(A)  complex  enhances  the  reaction  between  A  and  B  by  shifting  the  structure  of 
A  toward  what  it  is  expected  to  be  in  the  transition  state  for  the  neutral  reaction.  For  the 
above  reaction  that  is  toward  O  +  NO. 

Bond  Strength  vs.  Bond  Distance.  We  are  not  presently  able  to  estimate  how 
much  the  activation  energy  is  lowered  due  to  the  above  bond  lengthenings  and 
shortenings  but  we  can  estimate  how  much  the  N  -  0  bond  dissociation  energy  is 
reduced.  Figure  4  is  a  plot  of  the  logarithm  of  bond  strength  vs.  bond  distance  for  the 
N  -  O  bonds  in  several  molecules.  The  experimental  bond  lengths  in  NO,  N02,  and 
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N03  are  plotted  as  circles,  pluses,  and  squares,  respectively.  These  values  fall  on  a 
reasonably  straight  line  which  is  used  for  calibration.  The  calculated  values  for  the  N  - 
O  distances  in  Li+(N02)  and  N02  are  indicated  by  the  arrows.  The  estimated  difference 
in  the  bond  dissociation  energy  of  the  long  N  -  O  bond  in  Li+(N02)  and  the  bond 
dissociation  energy  of  N02  is  then  approximately  43  kcal  mol'1.  The  increase  in  bond 
strength  for  the  other  N  -  O  bond  in  Li+(N02)  is  about  the  same.  Therefore,  the  total 
bond  strength  of  N02  is  approximately  constant  in  Li+(N02)  and  N02.  Thus,  the 
increase  in  bond  length  translates  into  ca.  a  1/3  reduction  in  bond  strength.  This  must 
surely  translate  into  a  lowering  of  the  activation  barrier  as  well,  i.e.,  Ea'  <  Ea. 

G,  Components  of  the  Lithium  Cation  Affinities.  The  lithium  cation  affinity  can  be 
separated  into  four  main  attractive  components  -  electrostatic,  polarization,  charge 
transfer,  and  dispersion.26  For  first  row  bases  the  electrostatic  term  is  dominant;  for 
second  row  bases  both  the  electrostatic  and  polarization  terms  are  important.24  Since 
the  maximum  electrostatic  interaction  occurs  along  the  dipole  axis  of  a  molecule,27 
most  Li+  complexes  have  Li  •••  X  -  Y  angles  of  180°.102425  This  is  what  is  observed 
for  Li+(C02)  (Figure  3).  However,  ZLi  •••  O  -  N  in  Li+(NO)  is  not  180°,  although  the 
deviation  from  linearity  is  small  at  ca.  5°.  Moreover,  NO2  is  not  alligned  along  its 
symmetry  axis  in  Li+(N02)  (Figure  3).  The  latter  results  suggest  that  polarization 
effects  are  more  important  in  Li+(NO)  and  Li+(N02)  than  they  are  in  other  Li+ 
complexes.  In  fact,  the  magnitudes  of  many  of  the  parameters  used  to  measure  the 
size  of  the  polarization  contribution  are  larger  for  Li+(N02)  than  they  are  for  U+(C02). 
For  example,  forming  the  cluster  changes  the  O  -  X  (X  =  N,  C)  bond  lengths  more  in 
Li+(N02)  than  in  Li+(C02)  (Figure  3).  Furthermore,  the  increase  in  electron  density  on 
the  oxygen  bonded  to  lithium  and  the  decrease  in  electron  density  on  the  other  oxygen 
are  both  greater  for  N02  (0.182  e  and  0.186  e,  respectively)  than  they  are  for  C02 
(0.137  e  and  0.132  e,  respectively,  Figure  3)  when  the  complexes  are  formed.  In 
contrast,  the  electrostatic  contribution  to  the  lithium  cation  affinity  is  larger  tor  Li+(C02) 
than  it  is  for  Li+(N02).  The  Li  •••  O  distance  is  slightly  smaller  (Figure  3),  the  Is  orbital 
energy  of  the  Li+  accepting  oxygen  e1s(0)  is  less  stable  (-20.65782  au  vs.  -20.69138 
au),  and  the  charge  on  the  oxygen  in  the  isolated  base  is  Imore  negative  (Figure  3)  for 
C02.  The  compromise  between  the  two  components  accounts  for  the  similar 
magnitudes  of  the  electronic  lithium  cation  affinities  of  the  two  bases. 
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Reaction  Rate  Constant 

(cm3  s'1) 

Temp(K) 

Rate 

Enhancement0 

Ea 

kcal  mol' 

Li+(N02)  +  CO  ->  Li+  (C02)  +  NOa 

7(-12)b 

213-294 

% 

O 

CO 

28h 

Na+(N02)  +  CO  ->  Na+  (C02)  +  NO 

1  (“1 2) 

213-289 

1029  e 

28h 

Li+(N02)  +  H2  ->  products 

<1  (‘14) 

219 

- 

1 8h 

Li+(HBr)  +  0  -4  Li+  +  HBrO 

2(-10) 

178-219 

o 

Ol 

<0 

3.1h 

-4  Li+  +  Br  +  OH 

Li+(HBr)2  +  0  Li+  +  HBrO  +  HBr 

2(*10) 

178 

105  e 

3.1h 

Na+(HBr)  +  0  -4  Na+  +  HBrO 

i(-io) 

219 

2x1 O4  e 

3.1h 

->  Na+  +  Br  +  OH 

Li+(H2S)  +  0  ->  Li+(HSO)  +  Ha 

2(*10) 

219 

5x1 04  e 

4.3h 

-4  Li+  +  H2SO 

Li+(H2S)2  +  0-4  Li+(  H2SO  )  +  H2Sa 

2(-10) 

219 

5x1 04  e 

4.3h 

->  Li+  +  H2SO  +  H2S 

Na+(  H2S)  +  0-4  Na+(HSO)  +  Ha 

1(-10) 

178-219 

1.5x10s  e 

4.3h 

-4  Na+  +  H2SOh 

Na+(H2S)2  +  0-4  Na+(H2SO)  +  H2Sa 

2(*10) 

178 

3x1 05  e 

4.3h 

-+  Na+  +  H2SO  +  H2S 

Li+(N20)  +  O  -»  Li+  +  N2  +  02 

5(-13) 

219 

1025  0 

28h 

Li+(N205)  +  NO  -4  Li+(N204)  +  N02 

1.2(-11)d 

247-303 

>109  1 

22( 

Na+(N205)  +  NO  -4  Na+(N204)  +  N02 

3(-13)d 

303 

2:4x1 06  1 

22f 

K+(N205)  +  NO  -4  K+(N204)  +  no2 

^5(-14)d 

218 

- 

22* 

Li+(03)  +  NO  -4  U+  +  02  +  N02 

6(-1  1 )d 

220-280 

2x1 O4  9 

2.99 

Na+(03)  +  NO  -4  Na+  +  02  +  N02 

6.5(-11)d 

163-243 

2x1 05  9 

2.99 

K+(03)  +  NO  -4  K+  +  02  +  N02 

4.5(-12)d 

125-163 

2x1 05  9 

2.99 

aFor  detailed  information  on  the  products  see  text.  b7(-12)  means  7x1  O'1 2.  CThe  rate 
constant  in  the  presence  of  the  alkali  ion  divided  by  the  rate  constant  for  the  corresponding 
neutral  reaction  at  the  coldest  temperature  listed.  ^Data  from  Rowe  et  al.2  eNeutral  rate 
constant  calculated  from  the  tabulation  of  Kerr  and  Moss.20  *  A  limit  to  the  neutral  reaction 
rate  constant  can  be  obtained  from  the  thermal  decomposition  of  N205  in  the  presence  of  NO. 
The  data  of  Viggiano  et  al.18  are  used  here.  SNeutral  rate  constant  from  Birks  et  al.19 


Molecule  Point  Basis  HF  MP2 

group 


Li+ 


O 

NO 

C02 

N02 

LiO+ 


Li+(NO) 

Li+(N02) 

(II) 

U+(C02) 


6-31 G* 
6-31 +G* 

6-31 G* 
6-31 +G* 

6-31 G* 
6-31 +G* 

Dxh  6-31 G* 

6-31 +G* 

6-31 G* 
6-31 +G* 

6-31 G* 
6-31 +G* 

Cs  6-31 G* 

6-31 +G* 

Cs  6-31 G* 

6-31 +G* 

6-31 G* 
6-31 +G* 


-7.23554  a 
-7.23554  b 

-74.78393  a 
-74.78676 

-129.24788 

-129.25165 

•187.63418 

-187.63879 

-204.03149 

-204.03833 

-82.02031 

-82.02164 

-136.50479 

-136.50678 

-211.30064 

-211.30487 

-194.90516 

-194.90799 


-7.23554  b 
-74.88529 

-129.56539 

-188.11216 

-204.56557 

-82.12072 

-136.81482 

-21 1 .82729 


MP3 


-7.23554 

-74.89852 

-129.56347 

-188.09832 

-204.54433 

-82.13396 

-136.81500 

-211.80932 


-195.37855 


-195.36635 


Complex 

Basis 

HF 

MP2 

MP3 

LiO+ 

6-31 G* 

0.53 

6-31 +G* 

-0.41 b 

-0.069  b 

-0.063  b 

Li+(NO) 

6-31 G* 

13.41 

6-31 +G* 

12.29 

8.72 

10.03 

Li+(CO) 

6-31 1  +G(2d,2p) 

14.5  c 

Li+(N02) 

6-31 G* 

21.09 

(II) 

6-31 +G* 

19.45 

16.43 

18.48 

Li+(C02) 

6-31 G* 

22.24 

6-31 +G* 

21.12 

19.36 

20.39 

19.36 


20.39 


MICROWAVE 

DISCHARGE 


DIFFUSION  ROOTS  DIFFUSION 


It  would  be  useful  to  carry  out  theoretical  calculations  on  some  of  the 
Na+  clusters  and  more  of  the  Li+  clusters  studied  experimentally  in  this  work 
to  determine:  1)  additional  lithium  and  sodium  cation  affinities,  which  are  not 
known  experimentally,  2)  whether  the  structural  rearrangements  found  on 
complex  formation  are  a  universal  phenomenon,  and  3)  whether  the  relative 
rates  of  the  lithium  and  sodium  reactions  correlate  with  the  observed 
structural  rearrangements. 
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Collisional  quenching  of  NH  generated  from  N2/H2 
mixtures  by  high  energy  electron  impact  dissociation. 

Summary  Report. 


A  study  of  collisional  quenching  of  the  NH 

chemiluminescence  generated  by  the  high  energy  electron 
bombardment  of  mixtures  of  H.  and  excess  N.  has  been 
undertaken.  A  model  mechanism  has  oeen  formulates  for  the 
quenching  process.  The  pathway  of  imidogen  formation 
involves  several  steps:  primary  ion  formation  via  electron 
stripping  reactions 

N-  ♦  e”(40kev)  -->  N  +  ♦  e“(s) 

♦  e"(40kev)  -->  H*  ♦  e"  ( s ) 

generating  secondary  electrons  which  interact  more 
3trongly  with  the  initial  species  leading  to 

dissociation. 

N  ♦  e’  — >  2N(2S,2D,2P) 

+  e“  -->  2H(  S,  P) 

Of  the  dissociation  species,  only  N  D  has  the 
energetic  and  mechanistic  characteristics  necessary  to 
generate  NH(v<4). 

On  beam  termination  three  processes  dominate  the 
kinetics : 

production  of  v  ibr at ional ly  excited  NH 

N(2D)  ♦  H2  — >  NH ( v )  ♦  H  [1] 

2 

quenching  of  N  D  by  molecular  nitrogen 
N(2D)  ♦  N2  -->  N2«  ♦  N  [2] 

collisional  quenching  of  imidogen 

N H ( v )  ♦  N  — >NH  >  Np*  [3] 

NH ( v )  +  H2  — >NH  ♦  Hg*  [4] 

2 

If  the  assumptions  are  made  that  N  D  production  is  a 
consequence  of  electron  Irradiation  and  that  the  population 
of  excited  imidogen  is  proportional  to  that  of  NZD  during 
bombardment,  the  quenching  rates  of  imidogen  can  be 
deduced. 

The  rate  expression  for  imidogen  decay  for  this  model 
can  be  formulated  as: 

[NH(v)]/[NH(v) ]o  =  {-Aexpl-Bt)  ♦  Be xp { - At } } / { B-A } 


where  A  =  k1[H2]  ♦  k2  [  N2  ]  and  B  =  k3[N2]«-  k  4  C  H  ^  3  . 

The  concentrations  of  N2  and  H2  are  experimental 
parameters  and  the  rate  constants  kl  and  k2  have  been 
have  been  previously  de t erm  ined  [  1  ] ,  thus  the  quenching 
rates  k3  and  k4  may  be  obtain  through  curve  fitting  the 
experimentally  determined  ratio  [ NH ( v ) ] / [ NH ( v ) ] o . 


To  this  end  a  series  of  experiments  have  been  at  the 
Air  Force  Geophysics  Laboratory  conducted  using  the  "hot" 
LABCEDE  apparatus  which  has  been  described  in  detail 
elsewhere.  [2 ]  Continuous  flows  of  N2  and  H-  admixtures  at 
controlled  pressures  were  subjected  to  the  Doobardment  of 
a  pulse  collimated  37  kev  electron  beam  running  in  an 
Ilf  duty  cycle.  Fluorescence  was  monitored  by  means  of  an 
off  axis  Michelson  interferometer  with  a  LN-  cooled 
InSb  detector.  Interferogr ams  were  recorded  for 
mixtures  of  0.05,  0.10,  and  0.20  torr  H_  with  N-  at 
pressures  of  7.5,  10,  15,  20,  25,  30,  35,  and^UO  torr.  Fast 
Fourier  transformation  was  performed  on  each  of 
these  inter ferograms  and  the  resultant  spectra  obtained. 
The  desired  ratios  [ NH ( v ) ] / [ NH ( v ) ] o  was  then  determined  by 
integration  of  the  relevant  region  of  the  spectra. 

Curve  fitting  routines  have  been  devised  to  fit  the 
ratio  to  the  model  rate  equations.  Various  schemes  have  been 
tested  including  simultaneous  two  parameter  fits,  single 
parameter  fits  using  literature  values  for  the  appropriate 
rate  constants,  and  fits  with  a  variety  of  noise  treatments. 
Thus  far  only  whole  band  analysis  has  been  considered 
where  the  total  population  of  NH  independent  of 

vibrational  state  has  been  measured.  From  these 
preliminary  treatments  the  quenching  rate  constants  have 
been  estimated.  The  rate  of  quenching  of_JUi  by  N_,  k3, 

has  been  found  to  be  of  the  order  of  1x10“  cmr/aec  while 
the  quenching  rate  by  H_,  k4,  is  of  the  order  of  3x10“ 

cnr/sec.  It  must  be  stressed  that  these  finding  are 
preliminary  estimates  especially  for  k4  and  may  be  revised 
as  a  more  thorough  analysis  is  completed.  It  may  prove 
desirable  to  take  inter ferograms  at  additional  pressures 
to  improve  the  fit  of  k  4 .  A  better  accounting  of  the 
system  time  constants  would  also  be  highly  desirable 
perhaps  through  fits  to  an  RC  filtered  rate  equation  and 
state  dependent  analysis  remains  to  be  performed. 
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LINE  COUPLING  IN  THE  FUNDAMENTAL  Q  BRANCH  OF  CG2 

AT  6 67  CM-1 


MICHAEL  HOKE 


ABSTRACT 


Line  coupling  has  been  studied  in  the  fundamental  Q  branch  of  C02 
at  667  CM-1.  Transmittance  data  has  been  collected  using  the  high 
resolution  BOMEM  Fourier  transform  interferometer  at  the  National 
Bureau  of  Standards  Laboratory  in  Washington  D.C.  In  the 
simplified  analysis  line  coupling  coefficients  y  are  included  in 
the  line  shape  by  adding  a  term  for  each  line  proportional  to  y 
and  asymmetric  about  the  line  position.  The  coupling  coefficients 
are  given  in  terms  of  the  off  diagonal  elements  of  the  relaxation 
matrix  w(J,k).  These  elements  were  determined  from  the  collision 
broadened  half  widths  using  a  three  parameter  fitting  law 
expression.  The  agreement  between  the  observed  and  calculated  data 
at  the  pressures  for  which  data  was  recorded  (6*40,  380  and  190 
Torr)  was  good. 
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I)  INTRODUCTION  AND  OBJECTIVES 


Line  coupling  is  an  effect  which  is  increasingly  being 


recognized  as  important  in  the  interpretation  or  molecular 


spectra.  Also  referred  to  as  line  interference  or  rotational 


relaxation  the  effect  is  well  known  in  the  spectrum  o^  02  in  the 


microwave(l,2,3) •  The  effect  has  also  been  observed  in  Raman 


spectra  of  N2  and  C0(4)  and  N0(5)  and  more  recently  in  laser  diode 


absorption  spectra  of  C 02(6).  Line  coupling  has  been  discussed 


theoretically  in  the  literature ( 7 , 8 , 9 ,10) ,  and  there  has  been 


Increasing  interest  in  developing  computationally  efficient  and 


accurate  schemes  to  model  the  ef feet ( 1 , 4 , 6 , 11 , 12) .  We  have  applied 


a  simple  scheme,  after  Smith(3),  similar  to  Rosascod)  and 


Strow(6)  to  the  fundamental  Q  branch  of  C02.  This  0  branch  is  of 


particular  interest  because  the  lines  are  close  together  and  the 


coupling  large  but  also  because  of  its  use  in  remote  sensing  of 


the  atmosphere.  Recently  Armstrong(  13 )  and  Braund^)  have 


demonstrated  through  calculations  that  line  coupling  is  a 


significant  effect  in  this  0  branch.  What  follows  is  a  brief 


description  of  the  modeling  scheme  and  some  preliminary  results. 


II)  THEORY.  EXPERIMENTAL  DETAILS  AND  DATA  ANALYSIS 


The  expression  for  the  absorption  coefficient  K(v)  can  be 


written  formally  in  terms  of  the  spectral  density  function  as(3) 


k(v)  =  *Jtt2  n  {1  -  e"^hCy}  v  Im{I(v)} 


Eq.l 


with 


I(v)  =1  l  d . <  j|{(v-v°)  -  iPw } “ 1 | k  >d  p. 
*  jk  J  K 


Ea.2 


-  1..JW  ..O.  S!*.  a*  k  A  .  «  -  h  _  a* 


Here  <jjv|k>=v#I  where  v  is  the  frequency  and  I  the  identity 
matrix.  < j | v 0 | k>=v ( j ) *6 ( J , k )  a  diagonal  matrix  of  transition 
frequencies.  If  the  matrix  < j | w | k>=w( j ,k)  ,  referred  to  as  the 
rate  or  relaxation  matrix,  were  diagonal  the  impact  line  shape 
would  be  obtained.  The  Liouville  snace  states(8)  <j |  and  |k>  are 
products  of  the  upper  and  lower  Schrodinger  states  for  line  j  and 
line  k.  The  diagonal  elements  w(j,j)  of  the  rate  matrix  are  the 
impact ( Lorentzian)  half-widths.  The  off  diagonal  elements  w(j,k) 
connect  line  j  to  line  k.  Also,  in  Eq . 1  n  is  the  column  count 
[molecules/cm**2] ;  in  Eq . 2  the  d’s  are  reduced(3)  dipole  moment 
matrix  elements  and  p(k)  is  the  Boltzmann  factor  for  the  lower 
state  of  line  k. 

In  general  a  large  number  of  lines  are  usually  counled  and  to 
avoid  calculating  the  inverse  {v  -v0-iPw}-l  at  each  freouency  v 
Smith(3)  and  Rosenkranz ( 1 )  have  applied  a  perturbation  technique 
to  the  matrix  {v0+iPw}.  Writing  the  eigenvalues  and  eigenvectors 
as  power  series  with  pressure  as  the  expansion  parameter  Eq.3  is 
obtained,  for  the  Infrared,  to  first  order  in  pressure. 

Pwkk  P  (v-vk)  Yk 

- 4. - 

(v-vk)2  +  (Pwkk)2  (v-vk)2  +  (Pwkk)2  J  Ea . 3 


The  first  term  is  the  impact  line  shape.  In  the  second  term, 
asymmetric  with  respect  to  transition  frequency  v(k),  the 


quantities  y(k) 


I  (dj/dR)  wjk/(vk-vj) 
J  *  k 


Eq  . 


contain  the  contributions  from  line  coupling,  one  y  for  each  line. 


The  quantities  (  p(k)*d(k)*#2  )  and  (  d(j)/d(k)  )  can  be  related 
to  the  line  intensities  which  are  included  in  the  AFGL  Atmospheric 
Absorption  Line  Parameter  Compilation( 15) . 

The  elements  of  the  relaxation  matrix  w(j,k)  are  triven  as 
sums  of  matrix  elements  of  the  S  or  collision  matrix,  each 
averaged  over  all  impact  parameters  and  veloc  it  ies  (  2 , 3 , 1*0  •  Each 
element  contains  contributions  from  elastic  reorientat ional  and 
inelastic  rotational  transitions.  If  it  is  assumed  that  the 
inelastic  transition  contributions  dominate,  then  a  simple 
representation  of  the  rate  matrix  employing  "fitting  laws"(12) 
has  previously  been  used  to  model  molecular  spectra(  *4 , 5 ) .  We  have 
used  a  fitting  law  expression  (12): 


In  addition  if  only  inelastic  transition  rates  are  considered  the 


sum  over  the  rows  of  the  rate  matrix  of  a  given  column  is  zero, 


5  V  =  0 


Eq  .  8 


This  expression  allows  the  off  diagonal  elements  of  a  column  to  be 
related  to  the  diagonal  element,  which  is  the  Lorentz  width 


V  *  °k  *  - 


k  =  1  ,N 


Eo .  9 


Eq . 9  (with  Eqs.5,6)  is  used  to  estimate  Al,  A2  and  A3  which  are 
adjusted  until  the  sums  (k=l,...,N)  best  approximate  the  Lorentz 
widths,  in  the  sense  of  least  squares.  Since  the  model  w(J,k) 
Eq.6  is  non-linear  in  the  parameters  A2  and  A3  an  iterative 
fitting  routine  was  necessary. 


Ill)  RESULTS  AND  DISCUSSION 

Eq.9  has  been  fit  to  the  recently  measured  N2-broadened 
C02  widths  of  Arie  et  al(l6).  These  are  the  widths  which  will  be 
incorporated  into  the  new  1986  version  of  the  AFGL  Line  Parameter 
Compilation( 17 ) .  As  Fig.l  shows  the  difference  between  the 
one  atmosphere  observed  and  calculated  half-widths  is  less  than  5% 
over  most  of  the  range  of  M.  The  estimated  parameter  values  were : 
<A1>«0. 01211,  <A2>=0. 3272,  <A3>=1 .079.  With  Al ,  A2  and  A3 
determined  the  rate  matrix  w(j,k)  can  be  calculated  and  from  it 
the  y(k)'s,  which  are  shown  in  Fig. 2  and  listed  in  Table  1  for  the 
fundamental  Q  branch.  With  the  y’s  determined  the  spectrum  can  be 
calculated  in  a  straightforward  way  from  Eqs.1,3.  It  can  be  noted 
that  in  Eq . 3  each  asymmetric  term  varies  with  v  as  ±l/v  (depending 


on  the  sign  of  y)  far  from  the  line  center.  The  sum  of  these  terms 
can  grow  with  v  unless(7): 


N 


l  dkpkYk  =  0 


Eq.10 


which  will  be  true  only  if  the  rate  matrix  elements  satisfy 
detailed  balance,  Eq.7. 

One  calculation  for  the  fundamental  0  branch  is  shown  in  Fig. 3 
corresponding  to  the  specified  sample  conditions  listed  there.  The 
upper  curve  includes  line  coupling.  The  lower  curve  has  no 
coupling;  the  absorption  coefficient  is  the  superposition  of 
impact ( Lorentzian)  line  shapes.  Included  in  both  calculations  are 
the  contributions  from  the  many  other  (mostly  weak)  lines  in  the 
spectral  region  which  are  not  coupled  to  the  0  branch  lines  but  do 
contribute  to  the  observed  absorptance.  The  line  parameters  for 
these  lines  were  taken  from  the  1986  AFGL  Line  Parameter 
Compilation( 17 )  as  were  the  line  positions  and  strengths  for  the 
0  branch  lines. 

Also  included  in  Fig. 3  is  an  observed  spectrum  corresponding 
to  the  sample  conditions  listed  there.  It  was  obtained  with  the 
DOMEM  Fourier  transform  interferometer  at  the  National  Bureau  of 
Standards ( MBS )  Laboratory  in  Washington  D.C.  The  sample  was  an  NBS 
standard  reference  gas  mixture  of  0.96^0+0. 0005%  C02  in  dry 
nitrogen  contained  in  a  10  cm.  cell  with  KBr  windows.  The  optical 
path  was  evacuated  or  purged  with  nitrogen.  The  signal  was 
collected  with  a  liquid  nitrogen  cooled  Mercury  Cadmium  Telluride 
detector;  three  hundred  high  resolution  interf erograms  were 
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coadded  over  a  twelve  hour  period  before  transformation  to  yield  a 
signal  to  noise  ratio  of  approximately  65  to  1.  The  resolution  Is 
0.012  CM-1  using  a  Hamming  apodization  function.  Low  resolution 
(0.05  CM-1),  high  signal  to  noise  ratio  (1000  to  1)  empty  cell 
background  spectra  were  also  collected  before  and  after  the  high 
resolution  data.  They  were  in  very  good  agreement  and  were  used  to 
produce  the  Transmittance  spectr  m  shown  in  Fig. 3.  To  check  for 
instrumental  distortions  in  the  data  several  small  spectral 
regions  including  relatively  isolated  low  J  fundamental  p  and  R 
branch  lines  were  compared  with  calculations.  Over  a  fairly  large 
(=25  CM-1)  spectral  range  including  the  fundamental  Q  branch  the 
background  was  found  to  be  shifted  by  an  almost  constant  amount  of 
four  percent.  This  shift  is  included  in  the  calculations  in  Fig. 3 
Also,  a  constant  frequency  shift  of  approximately  0.015  CM-1  was 
eliminated  by  adjusting  the  positions  of  several  P  and  R  branch 
lines,  including  P 2,  P4 ,  R2  and  R^  in  Fig. 3.  Otherwise  the 
agreement  between  observed  and  calculated  P  and  R  branch  lines  was 
found  to  be  very  good. 

Considering  the  Q  branch,  the  difference  between  the  observed 
and  calculated  spectra  (obs.-cal.)  in  Fig. 3  shows  that  the 
agreement  is  good.  Some  difference  is  found  near  the  base  line 
where  the  absorption  coefficient  is  large  and  changing  rapidly 
with  frequency.  Since  the  elements  w(j,k)  are  given  in  terms  of 
the  collision  half-widths,  Eq.9,  they  would  be  expected  to  scale 
as  the  half-widths  do  with  pressure(see  also  Ben-Reuven( 8 ) ) .  Other 
room  temperature  spectra  were  collected  with  similar  experimental 


conditions  but  N2  broadening  pressures  of  38O  and  190  Torr.  The 
observed  and  calculated  data  were  in  good  agreement.  One 
atmosphere  air  broadened  line  coupling  coefficients  were 
determined  from  the  values  in  Table  1  by  assuming  02  broadening  to 
be  85£  of  N2  broadening,  based  on  the  02  broadened  C02  widths  of 
Arie  et  al(l6)  and  Bulanin  et  al(l8).  These  values  have  been 
included  in  the  1986  version  of  the  line  parameter  listing  and  are 
about  97?  of  those  in  Table  1. 

Strow  and  Gentry(6)  have  also  studied  line  coupling  in  C02, 
in  the  11102-00001  band  of  the  principal  isotope  located  at 
1932.47  CM-1.  These  authors  fit  the  N2  broadened  widths  of 
Yamaoto  et  al(19)  to  determine  Al,  A2,  A3  and  their  values  (and 
consequently  their  w(j,k)'s)  are  slightly  different  than  ours.  We 
also  tried  using  these  same  widths  but  found  me  widths  of  Arie  et 
al(l6)  gave  a  better  calculated  fit  to  the  observed  fundamental  P 
and  R  branch  lines,  including  those  shown  in  Fig. 3. 

Cousin  et  al  have  proposed  that  line  coupling  is  a 
significant  effect  in  the  v3  fundamental  R  branch  near  the  band 
head  and  have  employed  a  symmetrized  two  parameter  fitting  law 
model  including  only  the  exponential  gap  law.  For  comparison  we 
applied  their  model  to  the  widths  of  Arie  et  al(l6).  The  estimated 
Al,  A2  and  A3  (  and  w(J,k)  )  values  were  different  than  those 
given  above,  but  the  y's  calculated  from  them  were  very  close  to 
those  reported  in  Table  1. 

No  data  is  available  from  which  the  temperature  dependence 
of  the  rate  matrix  elements  and  y's  can  be  determined.  This  will 
be  a  topic  of  future  study. 
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TABLE  1 


NITROGEN 

BROADENED  LINE 

COUPLING 

COEFFICIENTS  y(M) 

M 

y  (M) 

M 

y  (M) 

2 

2.09 

44 

-0.0562 

4 

0.411 

46 

-0.0544 

6 

0.110 

48 

-0.0527 

8 

0.00737 

50 

-0.0512 

10 

-0.0387 

52 

-0.0494 

12 

-0.0602 

54 

-0.0480 

14 

-0.0709 

56 

-0.0465 

16 

-0.0763 

56 

-0.0452 

18 

-0.0781 

60 

-0.0439 

20 

-0.0786 

62 

-0.0427 

22 

-0.077 6 

64 

-0.0415 

24 

-0.0760 

66 

-0.0403 

26 

-0.0744 

68 

-0.0392 

28 

-0.0725 

70 

-0.0382 

30 

-0.0704 

72 

-0.0371 

32 

-0.0681 

74 

-0.0362 

34 

-0.0659 

76 

-0.0353 

36 

-0.0639 

78 

-0.0345 

38 

-0.0619 

80 

-0.0336 

40 

-0.0599 

82 

-0.0328 

42 

-0.0580 

84 

-0.0322 
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ed  using  the  model  described 


FIGURE  2:  One  atmosphere  line  couplir. 

for  N2  broadening  of  C02  ealcul 
<A1>=0. 01211 ,  <A2>=0.5272,  <A3> 
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for  the  Geopyhsics  Scholar  program  supported  by  SCEEE 
under  contract  ending  Aug.  31.  1986 

The  primary  focus  of  the  research  at  the  Air  Force  Geophysics 
Laboratory  has  been  on  improving  the  resolution  of  infrared  spectral  lines 
from  Fourier  transform  spectroscopy  <FTS). 

The  theory  is  based  on  the  minimum-negativity  constraint.  The 
theoretical  calculations  based  on  this  constraint  were  carried  out  earlier 
by  the  author^.  The  attempt  was  made  to  minimize  in  a  least-squares  sense 
the  negative  values  of  a  discrete  bandlimited  function  by  the  fitting  of  a 
discrete  function  formed  from  the  high  Fourier  frequencies  only.  This 
resulted  in  a  set  of  nonlinear  equations  in  the  coefficients  of  these  high 
Fourier  frequencies  which  have  to  be  solved  for  each  individual  set  of  data 
treated. 

The  work  at  AFGL  has  been  devoted  principally  to  finding  numerical 
solutions  to  these  nonlinear  equations  that  reduce  both  computer  memory 
requirements  and  c omp u ta t i ona 1  time.  The  only  known  means  of  solution  for 
these  equations  at  present  are  iterative.  The  Newton-Raphson  iterative 
procedure  may  be  used  to  obtain  a  rapidly  converging  solution  the  nonlinear 
equations.  However,  the  memory  of  the  computer  is  heavily  taxed  for  the 
large  data  sets  due  to  the  large  number  of  matrix  elements  involved. 

Other  matrix  have  proven  impractical  for  the  same  reason.  Incidentally, 
one  of  the  major  problems  with  the  maximum  entropy  methods  (MEM)  is  the 
large  number  of  matrix  elements  required  for  the  larger  data  sets. 

The  methods  of  solution  for  these  equations  that  have  proved  most  successful 
with  present  computers  are  the  various  relaxation  schemes  such  as  the  point 
successive  and  point  simultaneous  methods  of  successive  substitutions. 

The  point  successive  method  of  successive  substitution  was  adopted  by 
the  author  in  earlier  work*  as  the  equations  almost  always  converged  and 
converged  more  rapidly  than  the  point  simultaneous  method.  However,  the 
excessive  amount  of  cpu  time  needed  for  computing  the  larger  data  sets 
motivated  the  search  for  faster  methods. 

The  point  simultaneous  methods  are  being  explored  in  the  current 
research  being  done  at  AFCL  because  of  important  advantages  found  to  have 
been  possessed  by  them.  One  of  the  most  important  of  these  is  that  the 
fast  Fourier  transform  (FFT)  computational  algorithm  may  be  used  to 
calculate  the  summed  sinusoidal  products  in  the  equations  and.  of  course, 
to  perform  them  much  faster  than  may  be  computed  directly.  Another 
important  advantage  discovered  during  the  course  of  the  research  here 
was  that  there  are  ap p r o x l ma t i ons  that  may  be  employed  to  considerably 
reduce  the  overall  number  of  calculations  performed1. 

A  first  approximation  may  be  obtained  by  assuming  initial  values  of 
zero  for  all  the  unknowns  in  the  equations  except  the  one  whose  coefficient 
is  the  sum  over  the  squared  sinusoid  >  which  makes  this  coefficient  the 
largest  one  The  remaining  term  in  each  equation  is  easily  seen  to  be  a 
component  of  the  discrete  Fourier  transform  (DFT)  of  the  negative  values 
of  the  discrete  function  (  this  discrete  function  given  by  Fourier 
transforming  the  interferogram).  Thus,  the  a p p r o x ima t l ons  to  the  unknowns 
may  be  expressed  very  simply  by  factors  (which  are  not  constant,  but  vary 
with  each  unknown)  multiplying  the  FFT  of  the  negative  values  given  on  an 
iteration  A  reasonable  ap p r o x i ma t i on  to  the  interferogram  is  usually 
produced  by  this  operation  and  convergence  is  brought  about  for  many  data. 
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For  more  recalcitrant  data*  however*  better  approximations  are  needed. 
A  second  step  was  added  to  this  approx imation  by  substituting  the  unknowns 
approximated  by  the  above  procedure  back  into  the  equations  in  place  of 
those  unknowns  previously  assumed  to  be  zero.  This  brought  about  almost 
complete  convergence  for  all  data  considered.  However*  a  very  important 
discovery  made  in  the  course  of  this  work  was  that  not  all  of  the  terms 
are  needed  to  produce  convergence.  Only  a  subset  may  be  needed.  The 
subset  consisting  of  the  largest  values  of  the  terms  seemed  to  be  most 
effective  in  inducing  convergence.  Yet  total  cpu  time  was  not  reduced 
when  choosing  the  largest  values  directly.  This  is  because  the  decision¬ 
making  capabilities  of  the  computer*  such  as  expressed  by  FORTRAN  "IF" 
statements*  require  more  cpu  time  than  do  multiplications  and  other 
arithmetic  operations.  These  subsets  must  necessarily  be  chosen  by  less 
direct  means.  Much  time  and  effort  has  been  expended  in  this  research 
in  exploring  and  testing  the  various  ways  of  choosing  these  subsets  from 
the  total  number  of  terms  in  the  nonlinear  equations.  This  work  is  still 
ongoing. 

To  summarize  the  progress  at  this  point*  the  new  method  based  on  the 
point  simultaneous  solution  to  the  equations  reduces  computer  memory 
requirements  to  approximately  half  of  that  needed  by  the  earlier  method 
used  by  the  author*  and  reduces  cpu  time  to  about  one  fifth  of  that  of  the 
earlier  method,  for  reseanably  accurate  restorations  from  the  experimental 
data. 

For  future  development,  user  interaction  should  be  minimized  such 
that  the  average  technician  could  run  the  computer  programs  and  get 
reasonably  good  results.  At  present,  much  user  interaction  is  needed  to 
obtain  an  adequate  restoration.  Present  techniques  should  be  generalized 
and  further  developed  to  promote  the  smooth  processing  of  increasingly 
larger  data  sets*  for  both  one  and  two  dimensions.  Also,  various  other 
parameters  and  weights  for  the  subsets  and  other  terms  are  being 
adjusted  within  the  algorithm  in  efforts  to  find  the  combination  that 
results  in  the  most  complete  convergence  with  the  least  cpu  time. 

There  still  exists  a  problem  with  obtaining  complete  convergence  with 
the  experimental  data  when  using  the  approximations.  This  needs  to  be 
addressed  and  worked  out. 

The  data  employed  in  this  research  were  synthetic  infrared  OH  and 
COj  spectra  and  experimental  data  of  16384  points  in  the  near-infrared 
spectral  range  taken  with  a  rocketborne  f i e 1 d -w i d ened  interferometer 

1.  S.  J.  Howard*  J.  Opt.  Soc.  Am.  71,  819  <1981). 

2.  S.  J.  Howard,  Appl  Opt.  25*  1670  (1986) 
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Preface 


This  final  report  to  the  Air  Force  Geophysics  Scholar  Program  covers  the 
period  March  1,  1986  to  August  31,  1986.  It  consists  of  two  manuscripts  that 
will  be  published  in  geophysics  journals,  pending  clearance  by  AFGL  and 
review  by  the  journals  editors. 


ON  THE  TECTONIC  SIGNIFICANCE  OF  THE  "DELAYED  AFTERSHOCK"  SEQUENCE 
ASSOCIATED  WITH  THE  1977  SUMBA  -«RTHQUAKE 


Robert  McCaffrey 

Air  Force  Geophysics  Laboratory,  Hanscom  AFB,  Bedford  MA  01731 
and  Department  of  Earth,  Atmospheric,  and  Planetary  Sciences, 
MIT,  Cambridge,  MA  02139 


Abstract  The  "delayed  aftershock"  sequence  of  strike-slip  earthquakes 
following  the  great  Sumba  earthquake  of  August  1977  at  the  Java  trench  Is 
shown  to  have  occurred  within  the  overriding  plate  and  to  be  similar  in  depth 
and  source  mechanism  to  severs1  other  large  earthquakes  in  the  region  between 
1968  and  1983.  Historically,  strike-slip  earthquakes  are  the  most  energetic 
expression  of  plate  convergence  as  ur'Vrthrust  events  in  the  Sumba  region  are 
’-are  and  relatively  small.  The  October  1977  aftershock  events  are  consistent 
with  this  pattern  and  thus  a  direct  causal  relationship  between  the  mainshock 
and  these  strike  slip  aftershocks  is  not  indicated 


Introduction 


The  Sumba  earthquake  of  19  August  19  77,  with  a  seismic  moment  of 
approx imate ly  1031  Nm  and  a  normal  fault  mechanism,  indicates  that  some 
modification  of  the  subduct  ion  process  is  occurring  at  the  Java  trench 
^Cardwell  et  si  1981  Hanks,  1979),  but  debate  continues  on  whether  this 
earthquake  ruptured  only  the  upper  part  of  the  lithosphere  as  a  plate  bending 
event  Silver  and  ’crdon  1 98  1  .  Silver  et  a  1  1986  >  nr  ruptured  its  entire 

irkness  i  'liven  and  Kanamor  i  1 98d  l.vnnes  et  a,’  19°,S.  Spence.  1986. 

S’ewart,  1 9  '  R j  The  Angus*  earthquake  and  its  manv  aftershocks  were  followed 
t>v  a  "delayed  aftershock"  sequence  a  group  'f  earthquakes  that  occurred 
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approximately  200  ka  NW  of  and  starting  about  50  days  (in  early  October,  1977) 
after  the  mainshock  (with  a  combined  total  seismic  moment  of  approximately 
1019  Nra) .  In  contrast  to  the  earthquakes  beneath  the  trench,  that  had  E-W 
striking  normal- fault  mechanisms  and  N- trending  T  axes,  the  October  sequence 
had  strike-slip  fault  plane  solutions  with  N-trending  P  axes  (Dziewonski  et 
al.,  1981;  Fitch  et  al . ,  1981;  Spence,  1986).  Both  direct  (Spence,  1986)  and 
indirect  (Dziewonski  et  al.,  1981;  Fitch  et  al.,  1981)  relationships  between 
the  mainshock  and  the  "delayed  aftershocks"  have  been  inferred. 

The  October  "aftershock"  sequence  provides  a  rare  opportunity  to  study  the 
response  of  the  lithosphere  to  a  large  earthquake,  and,  in  this  respect,  a 
fundamental  question  is  whether  the  earthquakes  occurred  within  the  upper  or 
lower  plate.  Because  these  earthquakes  occurred  beneath  the  forearc  where  the 
cop  of  the  subducting  plate  is  still  less  than  70  km  deep  and  because  the 
hypocentral  depths  reported  by  the  International  Seismologlcal  Centre  (ISC) 
for  this  region  can  have  uncertainties  of  several  tens  of  kms  (McCaffrey, 
1986),  this  is  not  a  trivial  task.  In  a  recent  study,  Spence  (1986) 
constrained  focal  depths  by  picking  pP-P  times  and  relocated  the  events  by  the 
method  of  joint  hypocenter  determination.  Using  the  relocated  depths  and  a 
linear  extrapolation  of  the  top  of  the  subducting  plate  starting  at  the 
trench,  he  inferred  that  the  October  earthquakes  occurred  within  the 
subducting  plate. 

A  source  in  the  downgoing  plate  leads  to  the  problem  of  explaining  how  the 
maximum  deviatoric  stress,  originally  tensile  along  the  plate's  dip  direction, 
becomes  compress ional  in  the  same  direction.  The  explanation  offered  by  Spence 
(1986)  is  that  the  large  event  at  the  trench  caused  an  overshooting  of  the 
plate's  equilibrium  position  and  generated  a  strain  pulse  that  propagated  down 
the  slab.  If  the  Suaba  earthquake  did  indeed  relieve  all  tensile  stress  in  the 
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slab,  Chen  its  stress  drop  provides  an  upper  bound  on  the  magnitude  of  die 
ambient  deviatoric  stress  in  the  subducting  plate.  It  is  notable  in  this 
regard  that  several  workers  (Fitch  et  al .  ,  1981;  Hanks,  1979;  Silver  and 
Jordan,  1983)  have  remarked  that  the  Sumba  earthquake  may  have  had  one  of  the 
largest  stress  drops  calculated  for  a  large  earthquake,  on  the  order  of 
several  hundred  bars. 

The  purpose  of  this  paper  is  to  present  evidence  that  the  October  strike- 
slip  earthquake  series  occurred  within  the  overriding  plate  and  that  it  merely 
reflects  the  normal  response  of  the  upper  plate  to  convergence  in  this  region. 
The  fault  plane  solutions  and  depths  of  the  October  1977  earthquakes  were 
similar  to  those  of  several  large  earthquakes  in  the  region  between  1968  and 
1983  (Figure  1),  whose  combined  seismic  moment  was  approximately  3xlOl®  Nm.  In 
particular,  the  largest  of  the  October  1977  earthquakes  is  shown  here  to  be 
similar  in  mechanism  to,  but  smaller  in  seismic  moment  than  one  on  26  January 
1968,  an  event  not  directly  related  to  a  normal-fault  earthquake.  Likewise, 
the  other  strike-slip  earthquakes  were  not  associated  with  large  normal-fault 
earthquakes  at  the  Java  trench,  nor  did  they  occur  within  the  subducting 
plate.  While  it  is  possible  that  the  1977  Sumba  earthquake  to  some  degree 
triggered  the  strike-slip  earthquakes,  in  my  view  such  faulting  is  caused  by 
plate  convergence  and  would  occur  eventually  through  normal  strain 
accumulation.  Thus  a  direct  causal  relationship  between  the  normal-fault  event 
and  the  "delayed  aftershock"  zone  is  not  required.  More  importantly,  the 
scenario,  and  its  implications,  proposed  by  Spence  (1986),  in  which  the  great 
Sumba  earthquake  caused  the  previously  stretched  slab  to  go  into  downdip 
compression,  is  unsupported. 


Observation* 
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Earthquake  of  7  October  1977  (12hl0m)  This  earthquake  accounted  for  nearly 
half  of  the  seismic  energy  released  during  the  strike-slip  earthquakes  in 
October  1977  (Dziewonski  ec  al . ,  1981).  Vertical  component,  long-period  P 
waves  are  shown  in  Figure  2  along  with  synthetic  seismograms  for  centroid 
depths  of  12  and  22  km  below  sea  level.  The  source  depth,  strike,  dip,  rake, 
source  time  function,  and  seismic  moment  were  determined  by  least-squares 
inversion  of  the  amplitudes  of  the  observed  seismograms  (Nabelek,  1984; 
McCaffrey,  1986).  The  synthetic  seismograms  were  generated  using  a  source 
structure  constrained  by  a  nearby  seismic  refraction  profile  (profile  MSN  18 
(Curray  et  al.,  1977),  200  km  west  of  the  October  earthquakes)  (Table  1). 

The  solution  at  12  km  depth  (Figure  2a;  Table  2)  fits  the  data  much  better 
than  that  at  22  km  (Figure  2b)  both  visually  and  statistically;  the  difference 
in  variance  for  the  two  solutions  is  greater  than  a  factor  of  2.  Using 
waveforms,  Fitch  et  al.  (1981)  determined  depths  of  12  km  (below  the  sea 
floor)  both  for  this  event  and  for  the  second  largest  event  (16  October, 
21h09m)  in  the  sequence,  and  also  found  depths  of  10  and  11  km  for  two  other 
events  by  relative  relocations.  Adding  the  water  depth  (  =  4  km)  to  these 
values,  the  depths  below  sea  level  are  14  to  16  km. 

Depth  of  the  Plate  Interface  Thrust  earthquakes  with  shallow,  N-dipping 
nodal  planes  near  Sumba  Island,  150  kra  directly  to  the  east  of  the  October 
1977  earthquakes,  fall  between  24  and  30  km  depth  (Figure  1).  These  events 
mark  the  thrust  interface  between  the  Indian  Ocean  and  SE  Asian  plates. 
Chamalaun  et  al.  (1982)  infer  that  the  thickness  of  the  crust  north  of  Sumba 
is  approximately  24  km,  based  on  preliminary  analysis  of  a  seismic  refraction 
profile.  The  plate  interface  near  Sumba  must  be  at  least  this  deep.  Similarly, 
the  Moho  beneath  refraction  line  MSN- 18  (Table  1)  is  at  approximately  19  km 
depth  and,  because  sediments  in  the  Lombok  Basin  are  undeformed,  the  structure 
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revealed  by  MSN- 18  is  probably  that  of  the  upper  plate.  Therefore,  the  1977 


strike-slip  event,  being  well  above  this  interface,  was  certainly  in  the  upper 


plate.  I  also  point  out  that  all  of  the  events  relocated  by  Spence  (1986;  his 


Figure  3)  whose  largest  semimajor  axis  of  the  confidence  ellipse  is  less  than 


10  km,  are  in  the  depth  range  of  14  to  22  km;  i.e.,  also  within  the  upper 


plate . 


One  of  the  arguments  for  a  causal 


relationship  between  the  great  Sumba  earthquake  and  the  strike-slip 


aftershocks  is  the  apparent  historical  quiescence  of  the  region  in  which  the 


aftershocks  occurred  (Dziewonski  ec  al.,  1981;  Fitch  et  al.,  1981).  A  major 


point  of  this  paper,  however,  is  that  strike-slip  earthquakes  are  common  along 


the  eastern  Sunda  arc  and  are  explainable  in  terms  of  the  regions  tectonic 


development.  Elsewhere  (McCaffrey,  1986),  I  showed  that  strike-slip 


earthquakes  are  predominant  in  producing  shortening  within  the  forearc  and 


interpreted  this  as  a  direct  consequence  of  collision  with  the  Australian 


continent.  The  best  example  of  a  strike-slip  earthquake  within  the  upper  plate 


is  the  26  January  1968  (ra^-6.0)  earthquake,  one  of  the  5  largest  shallow 


events  (Mq>1019  Nm)  to  have  occurred  along  the  eastern  Sunda  arc  in  the  past 


25  years.  According  to  ISC  reports,  it  was  preceded  by  a  few  hours  by  one 


rn^-5.0  foreshock  and  was  followed  by  few  aftershocks  (only  7  earthquakes  are 


reported  in  the  following  year  within  100  km  of  this  event).  Fitch  (1972) 


determined  a  strike-slip  solution  for  this  earthquake  from  P  wave  first- 


motions.  Here  I  use  the  long-period  P  and  SH  waves  from  this  event  to 


constrain  its  source  depth  and  seismic  moment  and  to  refine  its  fault  plane 


solution . 


A  direct  compar ison  of  long -period  P  waves  for  the  1977  and  1968  earthquakes 


(Figure  3a)  reveals  great  similarity  in  waveforms,  implying  similarity  in 
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depths  and  source  mechanisms.  The  pulses  are  broader  for  the  1968  event,  most 
likely  because  of  a  greater  source  duration,  but  in  most  cases  they  correlate 
with  similar  pulses  in  the  1977  seismograms.  The  puises  for  the  1968  event  are 
delayed  in  time  with  respect  to  their  counterparts  in  the  1977  seismograms, 
indicating  a  deeper  source  and/or  a  longer  source  duration.  Note  that  several 
of  the  seismograms,  particularly  at  ANP,  MAT,  NA1 ,  and  CHG,  display  a  few 
cycles  of  high  amplitude  energy  after  the  first  cycle  of  the  P  wave.  These 
arrivals  probably  originate  from  source  structure  because  they  appear  at 
several  stations  and  in  both  events,  and  appear  to  have  the  same  apparent 
velocity  as  the  P  waves  Some  differences  in  the  seismograms  are  due  to  a 
difference  in  orientation  of  the  fault  plane  solution;  for  example,  SBA 
displays  a  dilatational  first-motion  for  the  1977  but  is  nodal  for  the  1968 
event  and  CHG  is  nodal  for  the  1977  but  compressional  for  the  1968  event 
Seismograms  and  the  best -fitting  source  mechanism  for  the  1968  earthquake 
(Table  2)  are  shown  in  Figure  3b;  this  solution  is  similar  to  that  determined 
by  Fitch  (1972).  The  1968  eartnquake  is  similar  in  centroid  depth  (16  *U  km) 
to  the  October  1977  event  (Table  2)  and  has  a  seismic  moment  il  2'  *0  31  xlQ19 
Nra)  that  is  greater  than  or  equal  to  that  of  the  entire  sequence  of  strike- 
slip  earthquakes  in  19  7  7  (Dziewonski  et  al  1 Q 8 1 

The  1968  earthquake  occurred  between  Sumba  and  Flores  Islands  very  close  to 
the  line  of  active  volcanoes  (.Figure  1)  Because  the  slab  here  Is  probablv 
nearly  100  km  deep,  there  can  be  no  doubt  that  this  cc"  ;  ,.iip  occur  red  in  the 
upper  plate  The  large  normal  fault  earth  ;  .jm1  at  the  'a-a  ♦  r  et.<  h  that  was 
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Strike-slip  faulting  within  the  upper  plate 


Other  Strike-slip  Earthquakes 
is  typical  of  the  stvle  of  deformation  for  this  part  of  the  Sunda  arc,  evident 
in  both  large  earthquakes  (Figure  1)  and  in  observed  llneations  in  morphology 
(Silver  er  al  ,  1983)  In  addition  to  the  1968  earthquake,  several  other  large 
shallow  earthquakes  between  1968  and  1983  (Figure  1  and  Table  3)  had  fault 
plane  solutions  similar  to  those  of  the  October  1977  sequence.  The  earthquake 
south  of  eastern  Sumbawa  (11  March  1982)  had  a  seismic  moment  similar  to  that 
of  the  largest  of  the  October  1977  events.  These  other  strike-slip  earthquakes 
were  not  associated  with  anomalously  large  slab  events.  Note  in  Table  3  that 
the  combined  seismic  moment  of  the  strike-slip  earthquakes  is  several  times 
that  of  the  underthrust  earthquakes. 

The  presence  of  strike -slip  earthquakes  beneath  the  forearc  north  of  the 
Timor  trough,  where  there  have  been  no  normal-fault  events,  indicates  a  lack 
of  spatial  correlation  between  normal-fault  earthquakes  at  the  trench  and 
s-rike-slip  earthquakes  in  the  forearc.  Thus  the  presence  of  a  normal-fault 
e •-  e  n  t  at  the  trench  is  not  a  necessary  condition  for  strike-slip  events  to 
o  c  c  u  r 


Interpretation  of  the  October  Earthquake  Sequence 

7V  e  <ibe  i  1  -'■> '  '  sequence  of  strike-slip  earthquakes  is  shown  here  to  have 
red  in  'he  upper  plate.  It  was  probably  triggered  by  the  large  normal- 
fault  earthquake  at  the  Java  trench,  nevertheless  the  earthquakes  are  a  direct 
'.resequence  of  plate  convergence  in  the  region 
I--,  figure  'he  geomerrv  of  the  plate  interface  Is  shown  and  superimposed  on 
: '  are  'he  positions  of  the  August  normal  -  fault  event  at  the  trench  and  the 
•r.ber  strike  slip  events  The  bathymetric  profile  is  taken  from  the  map  of 
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Hamilton  (1979)  along  118*E  and  the  slab  is  shown  dipping  northward  beneath 
the  forearc  at  an  angle  of  7*.  the  dip  angle  of  the  top  of  the  plate  observed 
at  the  trench.  At  this  constant  angle,  the  top  of  the  slab  reaches  24  km 
beneath  the  epicentral  region  of  the  October  earthquake  sequence.  This  depth 
probably  represents  a  minimum  for  the  plate  interface  as  discussed  above . 
Furthermore,  because  both  slab  geometry  and  lithospheric  bending  models 
indicate  that  the  dip  of  subducting  plates  increase  arcward,  the  linear 
extrapolation  used  in  Figure  4  will  also  cause  an  under-estimation  of  the 
depth  of  the  plate  interface.  Nevertheless,  even  this  minimum  depth  is  deeper 
than  the  source  depths  of  the  largest  of  the  October  1977  earthquakes  (Figure 
2)  and  of  the  other  large  aftershocks  (Fitch  et  al . ,  1981). 

There  is  little  doubt  that  the  Indian  Ocean  and  SE  Asian  plates  are 
converging  and  therefore  that  the  large  normal- fault  event  at  the  Java  trench 
represents  deformation  within  the  Indian  Ocean  plate.  If  slip  occurred  along  a 
45*  dipping  surface  (Given  and  Kanamori,  1980),  then  it  consisted  of  nearly 
equal  amounts  of  horizontal  and  vertical  slip.  However,  because  it  was  an 
intraplate  event,  it  is  unlikely  that  the  August  earthquake  accommodated  any 
net  motion  between  the  Indian  Ocean  plate  (south  of  the  Java  trench)  and  the 
SE  Asian  plate  (north  of  the  trench).  Thus  points  A  and  B  in  Figure  4  remained 
fixed  relative  to  each  other.  Furthermore,  as  Spence  (1986)  points  out,  there 
is  no  evidence  that  a  thrust  fault  was  activated  during  or  after  the  August 
earthquake  sequence,  an  observation  he  uses  to  infer  that  the  entire  Indian 
Ocean  plate  is  under  horizontal  tensile  deviatoric  stress.  This  lack  of  a 
* hrust  fault  also  means  that  the  interplate  thrust  zone  (Figure  4)  was 
if a--ttve  throughout  the  normal-fault  episode.  However,  because  the  normal 
‘  '  produced  divergence  between  points  A  and  C,  there  was  convergence 

..  -he  part  of  the  subducting  slab  downdip  from  the  normal  -  fault 
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earthquake*  (point  C)  and  tha  island  arc  (point  B)  that  suit  b*  taken  up  by 
strain  within  the  upper  plat*  The  sens*  of  this  strain  Is  the  sea*  as  that 
caused  by  long-tern  plate  convergence  in  the  region  (between  points  A  and  B 
chat  nomally  manifests  itself  by  strike-slip  faults 

Triggering  of  normal -fault  earthquakes  at  trenches  by  large  under thrust lng 
earthquakes  has  been  observed,  for  instance,  the  1965  Rat  Island  and  1960 
Chile  events  (Stauder,  1968;  1973)  It  is  reasonable  to  expect  reciprocity  In 
this  process;  that  normal -fault  earthquakes  at  the  trench  could  trigger 
underthrusting  events.  Because  the  selsmotectonlc  history  of  the  Sunda  arc 
shows  that  convergence  proceeds  predominantly  by  strike -slip  faults,  rather 
than  by  underthrusting ,  it  is  understandable  that  the  great  Sumba  earthquake 
would  cause  strike-slip  within  the  upper  plate  rather  than  lnterplat* 
thrusting.  The  Important  point  is  that  within  the  framework  of  our 
understanding  of  the  tectonics  of  the  Sunda  arc  collision  zone,  the  October 
1977  strike-slip  earthquakes  are  explained  in  a  simple  and  consistent  manner 
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TABLE  1  Observed  and  assuaed  velocity  ■tpictur*i 


Results  of  HSN  IS  at  10  2*S  115  )  *  E  -'urray  et  »1  19' 
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Structure  used  for  generation  of  synthetic  it  1  sinograms 
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TABLE  3.  Depths  and  seismic  moments  of  strike-slip  end  underthrusc  earthquakes 
from  the  eastern  Sunda  arc. 
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EARTHQUAKES ,  GRAVITY.  .AND  THE  ORIGIN  OF  THE  BALI  BASIN:  AN  EXAMPLE  OF  A 
NASCENT  CONTINENTAL  FOLD- AND -THRUST  BELT 


Robert  McCaffrey1  and  John  Nabelek2 
Department  of  Earth,  Atmospheric  and  Planetary  Sciences 
Massachusetts  Institute  of  Technology 
Cambridge,  MA  02139 


1  Nov  also  at  Air  Force  Geophysics  Laboratory,  Hanscom  AFB 
:  Nov  at  Lamont-Doherty  Geological  Observatory 

Abstrac t .  We  infer  from  the  bathymetry  and  gravity  field  and  from  t : 
source  mechanisms  and  depths  of  the  eight  largest  earthquakes  in  the 
region  that  the  Bali  Basin  is  a  downwarp  in  the  crust  of  the  Sundu  S1  t  : 
produced  and  maintained  by  shallow  thrusting  along  the  Flores  ha  i-  <  ■ 
zone.  Earthquake  source  mechanisms  and  focal  depths  are  in  f  »-• :  :  *  t  •  :  -  • 
inversion  of  long-period  P  and  SH  waves  and  from  short 
of  the  events.  Centroidal  depths  that  give  the  bes-  .•  • 

range  from  10  to  18  km  but  uncertainties  in  dep'r  i ; 
km.  The  P-wave  nodal  planes  that  dip  s.u.tr  «• 
parallel  to  the  volcanic  arc  and  are  ••  •  . 

Bali  Basin  beneath  it  The  p. > •  :  •  .  t 
crustal  features  inferred  ft  •-  ...  ; 

earthquakes  occur  in  the 
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zone.  Summation  of  seismic  moments  of  earthquakes  between  1960  and  1985 
suggests  a  minimum  rate  of  convergence  across  the  thrust  zone  of  4  ±2  mm/a. 

The  slip  direction  for  the  backarc  thrust  zone  inferred  from  the  orientation 
of  the  earthquake  slip  vectors  indicates  that  the  thrusting  in  the  Bali  Basin 
is  probably  part  of  the  overall  plate  convergence  as  it  roughly  coincides  with 
the  convergence  direction  between  the  Sunda  arc  and  the  Indian  Ocean  plate. 

The  presence  of  backarc  thrusting  suggests  that  some  coupling  between  the 
Indian  Ocean  plate  and  the  Sunda  arc  occurs  but  mechanisms  such  as  continental 
collision  or  a  shallow  subduction  of  the  Indian  Ocean  plate  probably  can  be 
ruled  out.  The  present  tectonic  setting  and  structure  of  the  Bali  Basin  is 
comparable  to  the  early  forelands  of  the  Andes  or  western  North  America  in 
that  a  fold-and- thrust  belt  is  forming  on  the  continental  side  of  an  arc- 
trench  system  at  which  oceanic  lithosphere  is  being  subducted.  The  Bali  Basin 
is  flanked  by  the  Tertiary  Java  Basin  to  the  west  and  the  oceanic  Flores  Basin 
to  the  east  and  thus  provides  an  actualistic  setting  for  the  development  of  a 
fold-and- thrust  belt  in  which  structure  and  timing  of  deformation  can  change 
significantly  along  strike  on  the  scale  a  few  hundred  kilometers. 

Introduction 

The  questions  of  how  and  under  what  conditions  fold-and- thrust  belts  develop 
in  the  forelands  of  continental  arcs  have  been  approached  largely  by  studies 
of  well  developed  examples  such  as  the  Andes  and  western  North  America.  In 
both  examples,  however,  the  along  strike  variations  in  the  style  of 
deformation  are  so  great  that  general  cause-and-effect  relationships  are 
difficult  to  infer.  Moreover,  in  only  one  of  the  cases  (the  Andes)  do  we  have 
first-hand  knowledge  of  the  dip  angle  of  the  subducting  plate  so  that 
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inferences  about  the  role  of  the  under thrusting  plate  in  driving  the 
deformation  in  both  examples  remain  equivocal.  Accordingly,  studies  of  other 
regions  where  backarc  thrusting  now  occurs  are  crucial. 

One  of  the  clearest  examples  of  present  day  backarc  thrusting  that  spans 
both  oceanic  and  continental  settings  is  the  long  thrust  zone  north  of  the 
Sunda  arc  between  Java  and  Wetar  [Hamilton,  1979;  McCaffrey  and  Nabelek,  1984, 
1986;  Silver  et  al.,  1983,  1986a, b;  Usna  et  al.,  1979]  (Fig.  1).  The  thrust 
zone  is  evident  in  two  segments:  the  Flores  thrust  zone  in  the  west  and  the 
Wetar  thrust  zone  in  the  east.  Both  dip  opposite  to  the  sense  of  subduction  of 
the  Indian  Ocean  -  Australia  plate  at  the  Java  Trench  and  Timor  Trough.  The 
Flores  thrust  zone  is  the  longer  and  more  developed  of  the  two  and  extends 
westward,  disappearing  beneath  the  Bali  Basin  [Silver  et  al.,  1983],  Both 
thrust  zones  have  gravity  anomalies,  accretionary  prisms  [Silver  et  al., 

1983],  and  thrust  earthquakes  (McCaffrey  and  Nabelek,  1984,  1986]  associated 
with  them.  So  far,  earthquakes  studied  that  are  clearly  associated  with  the 
backarc  thrust  zone  are  shallower  than  25  km  [McCaffrey,  1986]. 

In  this  paper  we  focus  on  the  western  section  of  the  Flores  thrust  zone  west 
of  Sumbawa  and  extending  into  the  Bali  Basin.  Here,  the  volcanic  arc  is  built 
on  the  thick  crust  of  the  southern  edge  of  the  Asian  margin  (Sunda  Shelf)  and 
thus  forms  a  continental  arc,  or  what  is  commonly  called  an  "Andean"  margin. 
The  Flores  thrust  zone  in  this  region  accommodates  the  thrusting  of  the  Sunda 
Shelf  in  the  backarc  beneath  the  volcanic  arc  and  thus  is  analogous  to 
intracontinental  thrust  zones  such  as  those  today  east  of  the  Andes  and  in 
western  North  America  in  the  Cenozoic.  To  the  east  the  Flores  and  Wetar 
backarc  thrust  zones  involve  oceanic  crust  of  the  Flores  Basin  and  Banda  Sea. 
We  examine  the  depths  and  mechanisms  of  the  eight  largest  earthquakes  from  the 
Bali  Basin  by  inverting  their  long-period  P  and  SH  waves  and,  in  some  cases, 
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short-period  P  waves.  In  addition,  we  compile  gravity  and  bathymetry  data  and 
present  maps  for  the  Bali  Basin.  We  then  discuss  the  constraints  these  data 
place  on  the  structure  and  evolution  of  the  Bali  Basin  and  the  causes  of 
backarc  thrusting,  and  discuss  the  implications  for  the  early  evolution  of 
fold- and- thrust  belts. 

The  Bali  Basin 


Setting  The  Bali  Basin  is  a  narrow,  easterly-elongated  basin  in  the 
southeastern  part  of  the  Sunda  Shelf,  with  water  depth  locally  exceeding  1.5 
km  (Figs.  1  and  2).  Its  southern  margin  is  formed  by  the  active  volcanic  arc. 
Along  strike  to  the  west  are  the  Tertiary  Madura  and  Java  sedimentary  basins 
and  to  the  east  is  the  deeper  Flores  Sea,  likely  underlain  by  oceanic  crust 
[Curray  et  al.,  1979;  McCaffrey  and  Nabelek,  1984], 


The  southeastern  Sunda  Shelf  is  covered  by  1  to  2  km,  and  locally  up  to  : 
km,  of  sediments  and  by  an  average  of  200  m  of  water  [Ben-Avraham,  1973;  Ben- 
Avraham  and  Emery,  1973].  It  locally  emerges  at  the  islands  of  Kangean  and 
Madura  in  an  E-W  trending  ridge.  During  most  of  the  Pleistocene  the  portion  of 
the  shelf  north  of  Java  was  a  landmass  that  submerged  only  recently  [van 
Bemmelen,  1949].  The  easterly  elongated  Java  Basin  received  sediments  from 
both  this  northern  landmass  and  the  volcanic  arc  to  the  south  [van  Bemmelen, 
1949;  Weeda,  1958].  The  eastern  Java  Basin  is  now  approximately  150  km  wide 
and  filled  with  up  to  6  km  of  Tertiary  sediments;  the  lowest  4  km  consist  of 
deep  water  facies,  while  the  upper  2  km  show  the  transition  to  the  present 
land  conditions  [Weeda,  1958].  We  know  of  no  direct  measure  of  the  thickness 
of  sediments  in  the  Bali  Basin,  but  we  suggest  that  there  may  be  approximately 
6  km,  based  on  gravity  data  and  the  shapes  of  the  seismograms  from  the 

98 


5 


earthquakes  (discussed  below) . 

Basement  rocks  of  the  eastern  Sunda  Shelf  range  in  ages  from  58  to  140  ma 
(K-Ar)  and  consist  of  terrigenous  and  volcaniclastic  metasediments,  and 
volcanic  and  granitic  rocks  ( Ben-Avrahara,  1973;  Hamilton,  1979],  The  basal 
sediment  layer  appears  to  be  of  Eocene  age.  The  sedimentary  cover  is  thin  or 
absent  over  basement  ridges  that  trend  either  NE  (the  Pulau  Laut  and  Meratus 
Ridges)  or  east  (the  Madura -Kangean  Ridge).  Ben-Avraham  and  Emery  [1973]  noted 
that  the  sedimentary  layers  on  both  sides  of  the  Pulau  Laut  and  Meratus  Ridges 
have  the  same  apparent  dip  as  the  ridge  flanks  and  do  not  thicken  away  from 
the  ridges.  This  geometry  indicates  that  the  ridges  formed  after  deposition  of 
the  strata  and  the  truncation  of  these  ridges  and  the  strata  at  the  seafloor 
suggests  that  uplift  has  ceased.  Conversely,  the  thickening  of  the  sedimentary 
layers  away  from  the  Madura -Kangean  Ridge  indicates  that  it  has  been  uplifted 
during  deposition  of  the  sediments,  and  its  present  exposure  at  the  two 
islands  suggests  that  uplift  may  be  continuing.  Ben-Avraham  and  Emery  [1973] 
interpreted  the  Pulau  Laut,  Meratus,  and  Madura -Range an  Ridges  as  the  roots  of 
anticlines  formed  by  NW-SE  and  N-S  compression. 

The  islands  of  Java,  Bali,  Lombok,  and  Surabawa  presumably  formed  by  the 
construction  of  volcanic  arcs  on  the  previously  passive  southern  margin  of  the 
Sunda  Shelf  [Hamilton,  1979],  Along  the  southern  coasts  of  the  islands  is  a 
pre-Miocene  volcanic  arc,  while  the  present  volcanoes  are  situated  closer  to 
the  northern  coasts.  Van  Bemraelen  [1949]  noted  a  NV-SE  alignment  of  volcanic 
features  on  Bali  and  a  tendency  for  younger  eruptions  to  be  more  basaltic. 
Between  Flores  and  Pantar,  the  volcanic  islands  are  built  on  oceanic  crust, 
the  volcanoes  lie  along  the  southern  coasts,  and  the  earlier  arc  is  not 
evident . 
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Crustal  Structure  North  of  Lombok,  the  Flores  thrust  zone  has  the 
appearance  of  a  typical  oceanic  subduction  zone  in  seismic  reflection  profiles 
(Line  63;  Fig.  3).  The  thick  sedimentary  section  south  of  km  30  is  broadlv 
folded  and  the  steeper  limbs  of  the  asymmetric  folds  are  truncated  along 
predominantly  S-dipping  thrust  faults.  The  increasing  dip  with  depth  of  planar 
reflectors  beneath  km  55  indicates  that  sedimentation  and  deformation  have 
occurred  concurrently.  North  of  the  Flores  thrust  zone  (km  0  to  20)  the 
reflectors  seem  to  arch  down  into  the  thrust  zone.  A  comparison  with  Line  67 
(Fig.  3)  shows  that  the  style  of  deformation  changes  westward  and  that  the 
Flores  thrust  zone  loses  its  surface  expression  north  of  Bali.  Nevertheless, 
the  broad  folding  continues  and  involves  the  entire  sedimentary  section, 
indicating  the  presence  of  a  decollement  surface  at  depth  beneath  the 
sediments  [Silver  et  al.,  1983].  The  gravity  field  is  not  appreciably  less 
negative,  indicating  that  the  amount  of  crustal  thickening  is  similar  beneath 
both  lines.  Note  again  the  arching  of  the  otherwise  undisturbed  reflectors 
from  the  north  into  the  deeper,  southern  part  of  the  basin. 

Hamilton  [1979]  suggested  that  the  crust  of  the  Bali  Basin  was  transitional 
in  thickness  between  oceanic  and  continental,  based  on  water  depth  and  the 
assumption  that  it  was  in  isostatic  equilibrium.  Ben-Avraham  and  Emery  [1973] 
inferred  that  oceanic  crust  was  present  beneath  the  Bali  and  Madura  Basins 
because  of  their  positions  on  line  with  the  Flores  Basin.  Based  on  the 
following  evidence,  we  infer  that  the  crystalline  crust  beneath  the  Bali  Basin 
is  the  same  thickness  as  beneath  the  Sunda  Shelf,  and  we  further  surmise  that 
the  two  are  genetically  identical.  First,  the  arching  cc  the  sediments  and  of 
the  seafloor  from  the  north  into  the  basin  (Figs.  3  and  U)  indicate  that  the 
basin  subsided  after  deposition  of  the  sediments.  The  continuation  of  the 
arching  into  the  deepest  part  of  the  basin  indicates  that  all  of  the 


bathymetric  relief  can  be  explained  by  downbowing  of  the  crust.  Second,  the 
minimum  in  the  gravity  field  over  the  Bali  Basin  is  displaced  30-40  km  south 
of  the  minimum  in  the  bathymetry,  indicating  that  isostatic  equilibrium  does 
not  hold  here.  More  likely,  the  basement  of  the  Sunda  Shelf  dips  to  the  south 
beneath  the  deformed  sediments  in  the  southern  Bali  Basin  and  the  deepening  of 
the  top  of  the  basement  is  not  accompanied  by  thinning  of  the  crust.  As  will 
be  shown,  the  magnitude  of  the  gravity  anomalies  is  consistent  with  the  idea 
that  the  gravity  and  bathymetric  lows  are  due  to  downbowing  of  the  crust  of 
the  Sunda  Shelf.  Finally,  the  bathymetric  profile  along  Line  22  of  Ben- 
Avraham  [1973]  (Line  B22;  Fig.  4)  displays  the  characteristic  shape  of  a 
flexed  lithospheric  plate.  The  estimated  flexural  rigidity  for  the  Sunda  Shelf 
of  5  to  10  xlO22  Nm  is  typical  of  some  oceanic  [Caldwell  et  al . ,  1976]  and 
continental  [Walcott,  1970]  regions  and  corresponds  to  an  effective  elastic 
thickness  of  18  to  23  km  (Fig.  4). 

Seismicity  Earthquakes  associated  with  subduction  of  the  Indian  Ocean  plate 
define  a  zone  that  dips  northward  from  the  Java  Trench  to  a  depth  greater  than 
600  km  beneath  the  Java  Sea  (Fig.  5).  We  note,  however,  that  in  this  region 
there  are  few  large  events  that  might  indicate  thrusting  of  the  Indian  Ocean 
plate  beneath  the  forearc,  whereas  several  normal  faulting  events  have 
occurred  beneath  the  Java  Trench,  including  che  -  8.1  earthquake  of  August 
1977  near  11*S,  118*E.  The  concentration  of  hypocenters  at  10*S,  117*E  is  due 
to  shallow  strike-slip  faulting  within  the  overriding  plate  during  October 
1977  [Fitch  et  al. ,  1981] . 

The  earthquakes  discussed  in  this  paper  occurted  at  shallow  depth  in  the 


overriding  plate  beneath  the  active  volcanic  arc  between  Java  and  Lombok,  and 
several  were  very  destructive.  Four  caused  damage  and  deaths  in  the  Bali 
region,  including  the  ra.  -6 . 1  event  of  *  luly  1976  (our  event  number  3)  that 


reportedly  killed  563  people.  In  the  sane  region,  earthquakes  in  1815  and  1917 
are  reported  to  have  killed  10,253  and  15,000  people,  respectively  [Ganse  and 
Nelson,  1981] . 

The  ambient  level  of  seismicity  in  the  epicentral  regions  of  the  large 
earthquakes  studied  here  is  very  low.  Since  all  events  from  the  ISC  bulletins 
are  plotted  in  Fig.  5a,  note  also  that  very  few  foreshocks  and  aftershocks 
with  magnitudes  greater  than  4.5  to  5.0  were  associated  with  these 
earthquakes . 

Analysis  of  Earthquakes 

Data  The  source  mechanisms  and  depths  of  the  earthquakes  listed  in  Table  1 
have  been  determined  by  matching  observed  long-period  P  (vertical  component) 
and  SH  waves  of  the  World-Wide  Standardized  Seismograph  Network  (WWSSN)  with 
synthetic  seismograms  (Fig.  6).  Each  earthquake  is  approximated  by  a  point 
source  having  a  double-couple  mechanism.  The  source  parameters  are  determined 
by  an  inversion  of  the  observed  data  in  a  least-squares  sense.  The  estimated 
parameters  are  the  strike  and  dip  of  one  P-wave  nodal  plane,  the  rake  angle  on 
the  plane  (using  the  convention  of  Aki  and  Richards  [1980]),  the  centroidal 
(average)  depth,  seismic  moment,  and  source  time  function  (i.e.,  source  time 
history).  We  will  refer  to  each  solution  by  the  strike,  dip,  and  rake  of  its 
N-dipping  P-wave  nodal  plane  because  its  strike  and  dip  angles  are  the  better- 
resolved  in  the  inversion.  This  plane  is  probably  the  auxiliary  plane.  The 
method  used  here  is  described  by  McCaffrey  [1986]  and  is  a  modification  of 
that  developed  by  Nabelek  [1984,  1985], 
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Only  seismograms  recorded  at  stations  at  epicentral  distances  between  30° 
and  90*  for  P  waves  and  30*  to  60*  for  SH  waves  are  used.  For  each  of  the 
earthquakes  presented  here,  between  12  and  29  seismograms  are  used  in  the 
inversion  (Table  1).  Data  are  the  amplitudes  of  the  observed  seismograms, 
digitized  at  0.25  s  intervals,  within  a  specified  time  window  starting  at  the 
direct  phase  (P  or  S)  and  extending  through  the  reflected  (pP,  sP  and  sS)  and 
other  converted  phases.  Attenuation  is  approximated  by  applying  a  causal 
attenuation  operator  (Futterman,  1962}  with  a  ratio  of  travel-time  to  average 
Q  of  Is  for  P  waves,  and  of  4s  for  SH  waves.  The  source  region  crustal 
structure  assumed  in  generating  the  seismograms  consists  of  a  water  layer, 
sediment  layer  (v^-3.2  km/s ,  v^-1 . 8  km/s ,  density  p-2300  kg/m3),  and  a 
halfspace  (v^-6.5  km/s,  vg-3 . 7  km/s,  p-2800  kg/m3).  The  P  wave  velocities  are 
based  on  nearby  seismic  refraction  results  reported  by  Ben-Avraham  and  Emery 
[1973],  Curray  et  al.  [1977]  and  Raitt  [1967],  The  thicknesses  of  the  water 
and  sediment  layer  (Table  1)  were  adjusted  to  match  specific  characteristics 
of  the  waveforms.  The  response  of  the  source  structure  was  calculated  by 
summing  a  finite  number  of  rays  [McCaffrey,  1986],  Structure  at  each  receiver 
was  assumed  to  be  a  half-space  (v^-6.0  km/s,  vs~3 . 4  km/s,  p-2700  kg/m3). 

Fault  Plane  Solutions  The  best-fitting  fault  plane  solutions  and  the 
observed  and  calculated  seismograms  are  shown  in  Fig.  6.  The  solutions  are 
also  plotted  in  both  map  and  cross-sectional  views  in  Fig.  7.  The 
uncertainties  in  the  source  parameters  given  in  Table  1  are  the  statistical 
uncertainties  and  underestimate  the  true  uncertainties.  An  evaluation  of  the 
actual  uncertainties  for  individual  events  is  presented  in  the  appendix  and 
the  general  results  are:  depth  ±5  km,  seismic  moment  ±50%,  strike  ±20*.  dip 
±5*.  and  rake  ±20*.  The  consistency  of  some  parameters  of  the  best-fitting 
solutions  among  the  events,  however,  indicates  that  the  uncertainties  might  be 
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less  than  those  just  cited. 

Events  1  and  2  occurred  near  the  NE  coast  of  Bali  in  1963.  Event  1  caused 
more  damage  on  Bali  although  it  was  deeper  and  smaller  in  seismic  moment  than 
event  2.  These  differences  imply  that  event  1  was  located  beneath  the  island 
while  event  2  was  offshore  (Fig.  7).  The  N-dipping  nodal  plane  of  event  2  is 
constrained  to  lie  between  stations  AAE,  SHI,  and  QUE  in  the  west  and  to  be 
near  GUA  in  the  east.  The  similarity  of  the  seismograms  (especially  from  QUE 
that  lies  very  close  to  the  nodal  plane)  for  the  two  events  suggests  very 
similar  mechanisms.  Of  the  events  studied  here,  these  have  the  shallowest  S- 
dipping  nodal  planes  (15*  and  14*). 

Events  3  and  4  occurred  near  the  north-central  coast  of  Bali  and  were 
separated  from  the  other  six.  Again,  the  N-dipping  nodal  plane  for  event  3  is 
constrained  by  the  northern  stations;  note  the  change  in  the  initial  P  pulse 
in  proceeding  azimuthally  from  CHG  (up)  to  HKC  and  ANP  (nodal)  to  SHK  and  HAT 
(up).  Event  4  is  small  and  occurred  in  the  surface  wave  coda  of  event  3.  The 
depths  of  these  two  events  estimated  from  the  short-period  P  seismograms 
(Table  1;  Fig.  8)  agree  to  within  2  km  with  the  depth  obtained  from  long- 
period  data.  This  agreement  is  compelling  evidence  that  the  long-period 
waveforms  do  have  sufficient  information  to  constrain  the  depths  of  these 
earthquakes  to  within  a  few  km. 

Events  5,  6,  7,  and  8  occurred  in  1979  beneath  the  strait  separating  Bali 
and  Lombok  Islands.  All  four  were  damaging  and  events  6,  7,  and  8  caused 
fatalities.  Note  the  strong  similarity  of  waveforms  among  the  events  with  the 
significant  exception  of  the  initial  half-cycle  at  stations  in  the  vicinity  of 
the  N-dipping  nodal  plane  on  the  focal  sphere. 

For  many  of  the  events,  long-period  seismograms  at  stations  near  the  N- 
dipping  P-wave  nodal  plane  display  a  very  sharp  initial  pulse  and  a  W-shaped 
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trough  (for  example,  see  NDI ,  CHG,  SHK,  and  MAT  for  event  3;  Figs.  6  and  9). 
This  shape  is  caused  by  the  presence  of  two  phases  with  negative  polarity 
(consistent  with  the  expected  polarity  of  pP;  sP  is  expected  to  be  positive) 
and  a  few  seconds  separation;  the  first  causes  the  sharp  turnaround  of  the  P 
and  the  second  produces  the  second  trough  of  the  V.  Note  in  Fig.  9  that  both 
the  long-period  and  short-period  records  at  HKC  are  nodal  (the  correlation  of 
later  phases  between  HKC  and  SHK  shows  that  HKC  has  no  direct  P  wave)  and  that 
the  first  pulse  observed  is  the  first  of  these  negative-polarity  phases.  The 
wave  shape  at  the  northern  stations  cannot  be  explained  by  adjusting  the 
source  depth,  the  source  time  function  (Fig.  9a  and  b) ,  or  the  thickness  of 
the  water  layer  within  limits  found  in  the  epicentral  region  (Fig.  9c).  A 
typical  Moho  interface  cannot  produce  the  observed  amplitudes. 

In  order  to  match  this  feature  of  the  observed  waveforms,  in  addition  to  a 
water  layer,  a  thick  layer  with  low  seismic  velocities  comparable  to  those  of 
sediments  was  required  in  the  source  region  (Fig.  9d,  e,  and  f ) .  In  this 
case,  the  first  negative  pulse  is  a  p-wave  reflection  from  the  bottom  of  the 
sediment  layer,  while  the  second  is  a  combination  of  the  p-wave  reflections 
from  the  tops  of  the  sediment  and  the  water  layer.  The  thickness  of  the 
sediment  layer  must  be  such  that  it  would  produce  the  time  difference  between 
the  phases,  while  the  impedance  contrast  at  its  base  must  be  large  enough  to 
produce  the  observed  amplitudes.  The  comparison  of  the  synthetic  seismograms 
produced  for  a  range  of  structures  (with  their  best-fitting  solutions) 
indicates  that  about  6  km  of  sediments  are  required  by  both  the  long-  and 
short-period  P  waves  (Fig.  9).  This  thickness  of  sediments  is  probably  only 
found  beneath  the  Bali  Basin  and  not  beneath  Bali  Island  or  Lombok  Strait, 
where  the  earthquakes  are  located  using  P-wave  arrival  times.  Thus  it  is 
likely  that  the  earthquakes  occurred  more  to  the  north  than  their  computed 
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locations  indicate.  Similar  mislocations  were  inferred  for  the  1978  earthquake 
on  the  Flores  thrust  zone  based  on  waveform  data  [McCaffrey  and  Nabelek,  1984. 
and  for  several  large  events  near  Timor  based  on  locally  recorded  arrival  time 
data  [McCaffrey  et  al . ,  1985]. 

Fault  plane  solutions  for  events  1-4  were  determined  from  P-  and  S-wave 
first  motion  polarity  readings  in  previous  studies  (event  2  by  Fitch  [1970], 
and  events  1,  3,  and  4  by  Kappel  [1980]).  The  best  fitting  strikes  of  all  new 
solutions  fall  within  10*  of  269*  (possibly  indicating  a  smaller  uncertainty 
in  strike  than  the  20*  estimated  in  the  Appendix)  whereas  those  of  the  earlier 
solutions  ranged  from  276  to  290*.  All  earlier  solutions  were  characterized  by 
rake  directions  of  90*  (i.e.,  pure  thrust),  whereas  the  revised  solutions 
indicate  a  systematic  deviation  from  pure  thrust  (all  best  fit  rake  angles  are 
within  9*  of  78*).  The  horizontal  projections  of  the  slip  vectors  (i.e., 
normal  to  the  strike  of  the  auxiliary  plane)  are  all  within  10“  of  N.  The 
convergence  direction  in  the  Bali  Basin  indicated  by  the  slip  vectors  is 
essentially  Identical  to  the  Indian  Ocean  -  Eurasia  convergence  direction 
determined  by  Cardwell  et  al .  [1981]  (0*)  and  is  similar  to  that  predicted  by 
Minster  and  Jordan  [1978]  (20*). 

Centroid  Depths  The  best-fitting  centroid  depths  determined  by  the  waveform 
analyses  range  from  10  to  18  km  but  the  uncertainties  in  depths  allow  a  range 
of  7  to  24  km  (Appendix) .  These  depths  are  generally  much  shallower  than 
those  reported  by  the  ISC,  including  those  based  on  "pP"  readings  (Table  1), 
and  further  Indicate  the  unreliability  of  the  ISC  depths  for  this  region 
[McCaffrey,  1986],  The  new  source  depths  are  unequivocal  evidence  that  these 
earthquakes  are  related  to  deformation  of  the  overriding  plate  and  not  to 
interplate  or  slab  activity.  The  long-period  data  do  not  contain  enough  high- 
frequency  information  to  distinguish  between  the  N-dipping  and  S-dipping  nodal 
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planes  as  the  fault  plane 

Stress  Drops  A  striking  feature  of  the  westernmost  events  (3  and  4)  is  that 
their  source  time  functions  are  very  short  when  compared  to  any  of  the  other 
events  of  similar  seismic  moment  and  that  the  moment  of  event  3  is  10  times 
larger  than  that  of  event  5  which  has  a  similar  duration  (Table  1).  In  a 
(single  asperity)  circular  crack  model,  the  stress  drop  is  proportional  to 
M  /t1  ;e.g  ,  Boatwright,  1980’.  where  t  is  the  source  duration  (Table  1). 
Because  M^/t3  is  roughly  an  order  of  magnitude  higher  for  events  3  and  4  than 
for  the  others,  it  is  likely  that  these  two  events  had  a  much  higher  stress 
drop.  It  is  unlikely  that  the  difference  is  due  to  error,  as  this  would 
require  a  factor  of  2  error  in  duration  or  a  factor  of  10  error  in  moment. 

The  higher  stress  drop  for  events  3  and  4  relative  to  those  of  the  eastern 
events  denotes  some  change  in  the  properties  of  the  fault  zone.  The  geometry 
of  the  Flores  thrust  zone  suggests  that  the  amount  of  slip  that  has  already 
occurred  decreases  westward.  It  is  thus  possible  that  events  3  and  4  occurred 
in  crust  that  is  less  fractured  and  perhaps  stronger  than  that  to  the  east. 

The  1978  Flores  and  1977  Timor  earthquakes  that  occurred  on  the  more  developed 
part  of  the  Flores  thrust  zone  and  on  the  Wetar  thrust  zone,  respectively,  are 
also  of  the  low  stress  drop  type  [McCaffrey  and  Nabelek,  1984,  1986], 

Estimate  of  the  Minimum  Slip  Rate  To  calculate  the  minimum  slip  rate  for 
the  Bali  Basin  we  use  the  relation  Iu-EMo/^A  normalized  by  the  time  period 
[Brune,  1968];  u  is  the  average  slip  for  one  event,  the  rigidity  (3.8 

xlO10  N/m2)  and  A  is  the  area  of  the  fault  surface  on  which  the  earthquakes 
occurred.  The  total  moment  221^  released  by  the  eight  earthquakes  was 
2  ±1  xlO19  Nm  (the  uncertainty  of  50%  is  based  on  the  analysis  in  the 
Appendix).  The  fault  length  along  strike  Is  estimated  at  150  km  (Fig.  7a)  and 
its  dovndip  length  at  40  km  (Fig.  7b),  for  an  area  of  6  xlO9  m2 .  The  time 
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period  is  25  years  (1960-1985).  From  the  relation  above,  the  total  slip  is 
0.09  ±0.04  m,  and  the  slip  rate  is  4  ±2  om/a  (rounded  to  the  nearest  mm/a) .  By 
a  global  fit  of  the  plates,  the  expected  convergence  rate  between  the  Indian 
Ocean  and  Asian  plates  near  Bali  is  approximately  70  mm/a  [Minster  and  Jordan, 
1978];  thus  the  rate  of  slip  across  the  Flores  thrust  zone  in  the  Bali  Basin, 
estimated  from  seismic  energy  release,  is  minor  compared  to  the  expected 
convergence  rate  between  the  major  plates.  We  emphasize,  however,  that  the 
lack  of  large  thrust  earthquakes  beneath  the  Sunda  forearc  [Kelleher  and 
McCann,  1976]  indicates  that  the  seismic  slip  rate  here  is  also  very  low 
compared  to  the  plate  motion  rate  and  note  that  the  concentration  of  thrust 
earthquakes  along  the  southern  margin  of  the  backarc  basin  is  also 
characteristic  of  the  Banda  arc  [McCaffrey  and  Nabelek,  1986] . 

Earthquake  Locations  The  eight  earthquakes  were  relocated  using  the 
P-tables  of  Herrin  [1968],  P-wave  arrival  times  from  ISC  Bulletins  and  the 
best-fitting  depths  determined  from  the  waveform  analysis  (Table  1).  The 
relocated  epicenters  are  shifted  northward  by  about  5  km  relative  to  the 
locations  determined  by  the  ISC  but  retain  their  original  pattern  of  grouping. 

Other  evidence  supports  the  relative  locations  determined  with  arrival  time 
data  and  gives  some  insight  into  the  absolute  locations  of  the  events.  Event  3 
produced  much  damage  in  western  Bali  and  was  felt  on  Java.  Its  location 
beneath  the  island  is  suggested  by  a  microearthquake  study  [Hayat  and 
Panjaitan,  1983]  in  which  shallow  seismic  activity  beneath  western  Bali  was 
recorded;  however,  the  level  of  activity  offshore  could  not  be  monitored.  The 
sequence  of  events  5  -  8  seems  to  have  migrated  from  the  vicinity  of  Lombok 
towards  Bali.  Events  5  and  6  caused  major  damage  on  Lombok  but  little  damage 
on  Bali.  Event  7  caused  considerable  damage  on  both  islands  whereas  the  damage 
caused  by  event  8  was  concentrated  on  Bali.  Event  8  produced  maximum  Modified 


Mercalli  Intensity  of  IX  due  east  of  Gunung  Agung  (the  southern  volcano  on 
Bali)  at  Culik  (8.3*S,  115. 6*E)  (Effendi  et  al . ,  1981],  suggesting  a  nearby 
epicenter.  This  constraint  on  the  epicenter  of  event  8  suggests  that  the 
epicentral  positions  shown  in  Fig.  7  be  shifted  approximately  10  km  northward. 
A  microearthquake  study  southeast  of  Gunung  Agung  [Wismaya,  1981]  also 
revealed  shallow  seismic  activity  beneath  Bali  Island.  Finally,  the 
requirement  of  a  thick  sediment  layer  in  the  source  structure  to  match  the 
waveforms  coupled  with  the  constraints  on  structure  from  the  gravity  data 
discussed  below,  suggests  that  the  earthquakes  were  farther  north,  beneath  the 
basin  where  sediments  are  more  likely  to  be  the  requisite  thickness. 

The  alignment  of  the  relocated  hypocenters  in  a  projection  onto  a  vertical 
plane  that  is  roughly  parallel  to  the  slip  direction  (Fig.  7b)  suggests  that 
the  S -dipping  plane  is  the  fault  plane.  Because  of  uncertainties  in  both  the 
epicenters  and  depths,  and  in  the  appropriate  projection,  this  independent 
determination  of  the  fault  plane  from  the  earthquake  positions  alone  must  be 
considered  somewhat  inconclusive.  Nevertheless,  the  sense  of  displacement  on 
the  Flores  thrust  zone  inferred  from  seismic  reflection  profiles  and  gravity 
data  suggest  a  S-dipping  plane,  and  the  fault  plane  solutions  suggest  that 
this  plane  dips  between  10*  and  30*.  No  significant  change  in  dip  angle  with 
depth  is  apparent. 
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Gravity  Interpretation 

Data  The  gravity  map  (Fig.  2)  and  Line  1  (Fig.  10a)  were  constructed  from 

gravity  data  compiled  by  C.O.  Bowin  (see  Bowin  et  al.  [1982])  and  the  gravity 

data  from  the  RAMA  12  cruise  in  1981  (Silver  et  al.,  1983].  All  gravity 

measurements  were  reduced  to  their  free -air  anomalies  with  respect  to  the  GRS 

67  formula  [e.g.,  Woollard,  1979].  The  gravity  values  on  Bali  Island  were 

reduced  to  simple  Bouguer  anomalies  using  the  infinite  sheet  approximation  and 

a  density  of  2700  kg/m3.  Gravity  and  bathymetry  values  were  averaged  over  5x5 

-4 

km  areas  and  contoured  at  20  mGal  (2  xlO  m/s2)  intervals  (Fig.  2). 

In  order  to  obtain  a  long  profile  over  the  Sunda  Shelf,  Bali  Basin  and  the 
Sunda  island  arc,  Line  1  was  constructed  by  collecting  all  measurements  within 
30  km  of  a  line  passing  southward  through  6'S,  116.41*E  at  an  azimuth  of  175' 
and  by  averaging  at  5-krn  spacings.  The  gravity  anomalies  along  Lines  63  and 
67  (Fig.  3)  were  taken  from  the  RAMA  12  data,  but  were  not  averaged.  In  the 

following  figures  all  profiles  are  referenced  to  the  same  x-axis  with  x-0  at 

6'S . 

Gravity  Map  (Fig.  2)  The  Bali  Basin  is  associated  with  a  low  in  the  free- 
air  gravity  field  that  follows  the  northern  coasts  of  the  volcanic  islands. 

The  low  reaches  below  -60  mGal  north  of  the  straits  between  Bali  and  Lombok 
and  between  Lombok  and  Sumbawa.  The  low  is  truncated  at  115°E  by  a  high  north 
of  western  Bali.  Over  the  Madura -Range an  Ridge  there  is  a  gravity  high  of  40 
to  60  mGal  that  is  likely  caused  by  the  combination  of  a  regional  high  and  the 
proximity  of  higher  density  basement  rocks  to  the  seafloor  [Ben-Avrahara  and 
Emery,  1973]. 

The  latitudinal  widths  of  both  the  gravity  and  bathymetric  lows  associated 
with  the  backarc  thrust  zone  appear  to  Increase  westward  from  Sumbawa  (Figs.  2 
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and  3).  We  note  Chat  throughout  the  length  of  the  Flores  thrust  zone  such 
variations  in  width  are  common  [Silver  at  al.,  1983 } ,  but  to  the  east  the 
regions  of  greater  width  are  usually  associated  with  more  thrusting  and 
folding  than  the  narrower  sections.  Beneath  the  Bali  Basin,  however,  the  zone 
of  deformation  is  wide  but  the  amount  of  shortening  in  the  sediments  appears 
to  be  small.  The  broadening  of  the  deformed  zone  probably  signals  a  change  in 
mechanical  behavior,  most  likely  in  the  sediments,  once  the  thrust  zone  steps 
onto  the  Sunda  Shelf. 

Isostaticallv  Balanced  Model  of  Line  1  All  gravity  and  bathymetric  lines 
show  a  strong  negative  in  the  gravity  anomaly  field  over  the  deep  parts  of  the 
Bali  Basin  and  then  an  increase  in  gravity  over  the  island  arc.  In  order  to 
determine  whether  or  not  the  crustal  structure  is  in  isostatic  equilibrium,  we 
generated  an  isostatically  balanced  model  (including  the  estimated  effect  of 
the  subducted  Indian  Ocean  plate)  for  Line  1  and  compared  the  calculated 
gravity  to  that  observed.  The  resulting  residual  (the  isostatic  anomaly)  is  a 
measure  of  the  amount  of  non-compensation  for  topographic  features.  The 
isostatically  balanced  crustal  model  and  its  gravity  effect  are  shown  in  Fig. 
10a. 

To  generate  the  balanced  model,  we  assumed  that  the  Sunda  Shelf  is  in 

isostatic  equilibrium.  This  is  a  safe  assumption  since  enormous  stresses  would 

be  required  to  keep  such  a  great  area  out  of  equilibrium.  For  the  shelf,  we 

assumed  a  crustal  thickness  of  25  km  (p^-2800  kg/m3).  The  average  water  depth 

of  0.1  km  (p^-1000  kg/m3)  and  sediment  thickness  of  1  km  (ps-2300  kg/m3)  are 

then  balanced  by  a  mantle  layer  thickness  of  1.9  km  (p^-3300  kg/m3)  using  the 

relation  h  (p  -p  )  +  h  (p  -p  )  +  h  (p  -p  )  -  0  where  h  is  the  thickness  of  the 
wwo  sso  mmo 

layer  whose  initial  is  in  the  subscript.  The  depth  of  compensation  (H)  is 
28  km  and  the  average  density  (pQ)  is  2800  kg/m3.  To  compensate  for  the 


changes  In  water  depth  the  crustal  thickness  h£  is  varied  locally  using  the 
relation  above,  and  the  relation  h  -H-h  -h  -h  .  The  sedimentary  layer  is 
included  in  order  to  employ  a  correct  density  for  many  of  the  bathymetric 
features  but,  because  its  thickness  is  kept  constant,  its  inclusion  has  a 
negligible  effect  on  the  long-wavelength  part  of  the  computed  anomaly. 

In  addition,  we  assumed  a  density  model  for  the  subducted  slab  of  the  Indian 
Ocean  plate  (Fig.  10a).  Part  of  the  motivation  to  include  a  slab  structure  was 
to  account  for  the  regional  gravity  high  over  the  Sunda  Shelf  [Ben-Avraham, 
1973],  The  position  and  shape  of  the  slab  are  based  on  earthquake  locations 
(Fig.  3).  A  density  contrast  with  respect  to  the  mantle  of  40  kg/ms  was  chosen 
to  match  roughly  the  amplitude  of  the  field  over  the  Sunda  Shelf  and  to  be 
consistent  with  the  observable  gravity  effect  of  slabs  in  other  parts  of  the 
world  [Grow  and  Bowin,  1975;  Molnar,  1977;  Watts  and  Taiwan!,  1975].  The 
gravity  effect  of  the  assumed  slab  structure  (shown  by  a  dashed  line  in  Fig. 
10a)  is  about  +35  mGal  over  the  northern  end  of  the  profile,  and  increases  to 
a  maximum  of  +70  mGal  over  the  island  arc.  Note  that  the  inclusion  of  the  slab 
as  shown  also  implies  isostatic  equilibrium  for  the  volcanic  islands;  they 
would  otherwise  have  a  strongly  negative  isostatic  anomaly. 

The  residual  gravity  of  »  -100  mGal  for  the  isostatically  balanced  model 
(isostatic  anomaly)  indicates  that  the  low  free-air  anomaly  observed  over  the 
Bali  Basin  is  not  an  effect  of  the  juxtaposition  of  crustal  elements  of 
different  thicknesses  (i.e.,  an  edge  effect).  We  conclude  that  the  crust 
beneath  the  Bali  Basin  does  not  thin  to  accommodate  the  increase  in  water 
depth;  instead  it  appears  to  be  deflected  downward,  and  the  void  is  filled  by 


water  and  sediments. 
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In  order  Co  estimate  how  much  the  crust 


beneath  the  Bali  Basin  is  downvarped,  an  isostatically  balanced,  constant 
crustal  thickness  model  of  the  Sunda  Shelf  was  perturbed  by  varying  the 
thickness  of  the  sediment  layer  until  agreement  with  the  gravity  data  was 
obtained  (Fig.  10b).  The  gravity  profile  for  Line  67  (Fig.  3b)  rather  than 
Line  1  is  used  because  this  line  is  the  closest  to  where  the  thrust 
earthquakes  occurred.  Again  the  densities  and  thicknesses  noted  above  for  the 
isostatically  balanced  Sunda  Shelf  were  assumed  for  this  profile  and  the  slab 
structure  was  included. 

To  match  the  gravity  data  the  following  perturbations  were  necessary 
(Fig.  10b):  (1)  thickening  the  sedimentary  layer  to  4  b  beneath  km  240,  (2)  a 
sharp  offset  in  basement  at  km  240,  (3)  a  small  ridge  in  basement  beneath  km 
200,  and  (4)  thinning  the  sediment  layer  beneath  the  southern  end  of  the  line. 
Thickening  of  the  sediments  beneath  the  Bali  Basin  was  necessary  to  produce 
the  large  gravity  negative.  While  the  amount  of  sediment  thickening  depends  on 
the  assumed  densities,  the  important  point  is  that  thickening  by  downbowing  of 
the  crust  is  required.  If  the  crystalline  crust  was  thinner  beneath  the  basin, 
then  even  more  downbowing  is  indicated.  For  the  densities  used  here  the  dip 
angle  of  the  basement  between  km  210  and  240  is  about  3* .  The  dip  angle  of  the 
major  thrust  surface  of  the  Flores  thrust  zone  between  km  20  and  35  In  Fig.  3a 
is  slightly  larger,  6*.  probably  because,  as  can  be  seen  in  Fig.  3a,  the 
thrust  plane  splays  from  a  probable  decollement  surface  at  depth  up  through 
the  sedimentary  section. 

The  sharp  crustal  offset  at  km  240  and  the  thinning  of  the  sediments  south 
of  that  point  was  demanded  by  the  sharp  gradient  in  gravity  from  km  240  to 
250.  Even  with  a  nearly  vertical  offset  in  basement,  the  calculated  gradient 
is  less  than  that  observed.  While  the  observed  gradient  may  indicate  a  greater 


densicy  contrast  between  the  crust  and  sediments,  the  lack  of  correlation 


between  the  gravity  field  and  bathymetry  In  the  strait  between  Bali  and  Lombok 
(Fig.  2)  suggests  that  It  Is  more  likely  due  to  structures  to  the  sides  of  the 
profile.  Note  that  the  contour  lines  of  bathymetry  are  nearly  tangential  to 
the  profile  line  at  its  southern  end  (Fig.  2)  and  thus  the  assumption  of 
structural  two-dimensionality  is  violated.  Finally,  since  the  crustal 
thickness  is  kept  constant,  deflection  of  the  Moho  is  implied,  but  we  point 
out  that  the  gravity  data  are  insensitive  to  such  details  of  the  shape  of  the 
Moho . 


Discussion 

Causes  of  Thrusting  In  the  section  of  the  Sunda  arc  from  Java  to  Wetar,  the 
convergent  margin  changes  from  subduction  of  old  oceanic  crust  (Indian  Ocean) 
beneath  a  continental  margin  (Sunda  Shelf)  in  the  west,  to  subduction  of 
continental  crust  (Australia)  beneath  an  island  arc  built  on  oceanic  crust 
(Banda  Arc)  in  the  east.  In  the  east,  most  of  the  deformation  of  the  island 
arc  can  be  attributed  to  collision  with  the  Australian  continent.  In  the 
Flores  section,  the  deformation  of  the  arc  is  probably  caused  by  interaction 
between  the  subducting  Australian  margin  and  a  thick  crustal  block  (Sumba) 
within  the  forearc.  At  the  western  end,  near  Bali,  however,  there  is  no 
obvious  local  cause  of  the  thrusting  north  of  the  volcanic  arc. 

One  explanation  for  the  presence  of  thrusting  in  the  Bali  Basin  offered  by 
Silver  et  al.  (1983)  is  lateral  propagation  of  the  Flores  thrust  zone  from  the 
east.  This  mechanism  requires  that  the  arc  structure  act  as  a  beam  when  it  is 
indented.  Because  the  arc  itself  is  cut  extensively  by  strike-slip  faults  and 
because  the  northward  displacement  of  the  volcanic  islands  north  of  Timor  and 


Sumba  appears  to  occur  only  locally,  we  feel  chat  lateral  propagation  of  the 
backarc  thrusting  Is  Inefficient  and  does  not  offer  a  viable  explanation. 
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Hamilton  [1979]  suggested  a  mechanism  for  the  thrusting  that  centered  on  the 
Intrusion  of  magmatic  material  and  required  no  external  forces  or  net 
shortening  across  the  arc.  He  observed  that  on  Java  "compresslonal  structures 
In  Neogene  materials  tend  to  arc  concentrically  around,  and  to  be  directed 
outward  from,  large  volcanic  edifices".  He  suggests  that  the  downbowing  that 
allowed  the  formation  of  sedimentary  basins  was  due  to  loading  of  the  crust  by 
volcanoes  and  that  subsequent  deformation  of  the  sediments  was  caused  "partly 
by  gravitational  flow  of  the  upper  crust  and  its  sedimentary  cover  in  response 
to  magmatic  loading,  and  partly  by  spreading  of  the  crust  in  the  magmatic  belt 
by  intrusions  within  it". 

The  distinction  between  the  mechanism  advocated  by  Hamilton  in  which  no  net 
shortening  across  the  volcanic  arc  structure  occurs  and  one  in  which 
shortening  does  occur  is  basic  to  our  understanding  of  island  arcs, 
collisions,  and  whether  or  not  island  arcs  can  reverse  subduction  direction. 
While  other  circumstances  may  apply  to  arcs  in  general,  we  favor  the  presence 
of  shortening  across  the  Sunda  island  arc  based  on  our  studies  here.  First, 
the  lower  crustal  depths  and  shallowly  dipping  nodal  planes  for  the 
earthquakes  beneath  the  Bali  Basin  are  evidence  that  basement  is  actively 
involved  in  the  thrust  faulting.  Second,  the  similarity  of  the  fault  plane 
solutions  despite  the  different  locations  of  the  earthquakes  relative  to  the 
active  volcanoes  suggests  that  the  thrusting  direction  is  not  controlled  by 
the  positions  of  the  volcanoes  but  is  due  to  more  uniform  slip  of  the  island 
arc  over  the  backarc  basin.  Finally,  the  most  critical  evidence  is  the 
locations  of  these  earthquakes  and  their  relationship  to  the  volcanoes  on  Bali 
and  Lombok  Islands.  The  locations  based  on  arrival  time  data  (Fig.  7)  and  the 
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intensity  reports  suggest  that  the  thrust  earthquakes  and  thus  the  thrust 
fault  zone  extend  beneath  the  active  volcanoes.  If  so,  then  mechanisms  for  the 
thrusting  that  rely  on  swelling  of  the  arc  by  magmatic  intrusion  from  beneath 
cannot  apply  and  regional  compression  of  the  volcanic  arc  must  be  invoked. 
Silver  et  al.  [1983]  also  rejected  the  magmatic  intrusion  mechanism  because 
the  backarc  thrusting  is  observed  where  there  are  no  volcanoes  and  because  of 
the  absence  of  obvious  extensional  features  within  the  arc  in  places  where 
backarc  thrusting  is  well -developed  and  significant  slip  (i.e.,  tens  of  kms) 
has  occurred  (e.g.,  Flores).  Fault  plane  solutions  for  shallow  earthquakes 
within  the  arc  and  forearc  are  predominantly  strike-slip  and  indicate  that 
they  are  in  fact  being  shortened  in  the  N-S  direction  [McCaffrey,  1986], 

What  causes  the  N-S  compression  across  the  arc?  The  roughly  parallel 
convergence  directions  between  the  Sunda  Shelf  and  the  island  arc  in  the  Bali 
Basin  indicated  by  slip  vectors  and  between  the  Indian  Ocean  plate  and  Asia 
[Cardwell  et  al.,  1981;  Minster  and  Jordan,  1978]  suggest  that  plate  motions 
are  driving  the  backarc  thrusting;  probably  through  coupling  of  the  plates 
beneath  the  forearc.  Such  coupling  could  be  produced  by  shallow  subduction  of 
the  Indian  Ocean  plate,  but  we  cannot  document  this  because  the  uncertainties 
in  the  earthquake  hypocenters  do  not  permit  resolution  of  the  subduction  angle 
of  the  Indian  Ocean  plate  at  shallow  depths  beneath  the  forearc  south  of  Bali 
and  Lombok.  There  are  other  indications,  however,  that  subduction  may  be  steep 
here.  The  subducting  oceanic  crust  is  very  old;  both  magnetic  anomaly 
lineations  (Larson,  1975]  and  upper  Jurassic  basal  sediments  in  Deep  Sea 
Drilling  Project  Hole  261  (at  1 3 * S ,  118'E)  indicate  an  age  of  about  late 
Jurassic  to  early  Cretaceous.  The  Lombok  Basin  is  the  deepest  section  of  the 
forearc  basin  (with  the  exception  of  the  Weber  Deep  far  to  the  east)  and  the 
Java  Trench  appears  to  be  slightly  deeper  south  of  the  Lombok  Basin  (Fig.  1). 


There  have  been  no  large  thrust  earthquakes  here  [Kelleher  and  McCann,  1976], 
but  there  seems  to  be  a  tendency  for  normal  faulting  events  at  the  trench 
(possibly  related  to  sharp  downbending  of  the  oceanic  plate)  [Cardwell  and 
Isacks,  1978],  including  the  Mg-8 . 1  Sumba  earthquake  of  1977.  The  foregoing 
observations  might  argue  for  a  steep  subduction  angle  and  very  weakly  coupled 
subduction.  In  fact,  considerations  of  such  factors  as  slab  dip  (at  greater 
depths  beneath  the  arc),  convergence  rate,  lithospheric  age,  and  maximum 
earthquake  size,  place  the  Java  trench  globally  among  subduction  zones  at 
which  backarc  spreading,  and  not  backarc  compression,  occurs  [Ruff  and 
Kanamori,  1980;  Uyeda  and  Kanamori,  1979]. 

Alternatively,  as  Silver  et  al.  [1983]  suggest,  and  we  favor  this  view,  the 
thrusting  in  the  Bali  Basin  may  be  due  to  collision  with  the  Roo  Rise,  a 
bathymetric  high  that  intersects  the  Java  trench  between  110'  and  115°E  (Fig. 
1)  and  probably  extends  to  the  NE  beneath  the  forearc  south  of  Bali.  Kelleher 
and  McCann  [1976]  inferred  that  the  Roo  Rise  resists  subduction  because  of  its 
buoyancy  and  is  responsible  for  the  lack  of  great  thrust  earthquakes  in  the 
forearc  south  of  Java.  Details  of  the  interaction  of  the  Roo  Rise  with  the 
forearc  are  lacking,  but  lateral  changes  in  forearc  structure  appear  to 
coincide  with  the  Roo  Rise  (the  forearc  outer  ridge  is  more  elevated  and  the 
forearc  basin  is  narrower  here  than  to  the  east  and  the  Java  Trench  appears  to 
be  slightly  indented  [Hamilton,  1979]),  indicating  that  some  modification  of 
the  subduction  process  is  taking  place. 

Analogy  with  the  Andes  and  western  North  America  In  the  forelands  of  both 
the  Andes  and  western  North  America  two  styles  of  deformation  are  dominant; 

(1)  the  thin-skinned  style  of  the  Canadian  Rockies  [Bally  et  al.,  1966;  Price, 
1981;  Price  and  Mountjoy,  1970],  Sevier  Desert  [Armstrong,  1968],  Bolivia,  and 
Argentina  [Jordan  et  al.,  1983]  and  (2)  block  faulting,  or  thick-skinned, 


style  of  the  Laramide  (Berg,  1962;  Prucha  et  al.,  1965;  Smithson  et  al . ,  1979] 
and  Pampeanas  [Jordan  et  al.,  1983]  ranges.  In  the  Andes,  while  the  thin- 
skinned  deformation  occurs  over  both  nearly  horizontal  and  shallow  dipping 
(=  30*)  subducting  slabs,  the  thick-skinned  deformation  is  limited  to  regions 
where  the  slabs  are  nearly  horizontal.  In  the  western  U.S.  the  development  of 
the  Laramide  basement  faults  was  accompanied  by  an  eastward  migration  of  the 
volcanic  arc  and  this  has  led  to  the  popular  view  that  the  slab  was  very 
shallowly  dipping  during  the  episode  of  Laramide  deformation. 

The  style  of  deformation  in  the  Bali  Basin  appears  to  be  intermediate 
between  that  of  the  thin-skinned  style  east  of  the  Andes  and  Rockies  and  the 
basement  block- faulting  of  the  Laramide  or  the  Pampeanas.  While  the 
earthquakes  examined  in  this  report  certainly  occur  within  basement,  and  not 
along  a  decollement  at  the  base  of  the  sedimentary  sequence,  they  are  also  of 
fairly  low  angle  and  the  associated  basement  faults  do  not  seem  to  extend  into 
the  sedimentary  layer.  The  sediments  respond  to  basement  shortening  by 
folding  and  are  probably  decoupled  from  basement  along  an  aseismic  shallow 
dipping  plane  that  steepens  into  basement  and  becomes  seismic  only  beneath  the 
southern  edge  of  the  Bali  Basin.  The  Madura -Kangean,  Meratus,  and  Pulau  Laut 
Ridges  are  basement  uplifts  that  sit  in  the  same  tectonic  position  as  the 
Laramide  and  Pampeanas  basement  uplifts  and  the  available  information 
indicates  that  they  are  the  roots  of  anticlines  formed  by  compression  [Ben- 
Avraham  and  Emery,  1973). 

The  Java-Bali  section  of  the  Sunda  Arc  is  a  likely  analog  for  the  initial 
stages  of  thrusting  on  the  east  side  of  the  Andes  (Jordan  et  al.,  1983]. 

Suarez  et  al.  (1983)  propose  an  evolutionary  scenario  for  the  Andes  that 
begins  with  the  Brazilian  Shield  thrusting  westward  beneath  the  volcanic 
chain,  forming  a  small  foreland  basin  and  uplifting  the  mountain  belt.  As  the 
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thrusting  and  uplift  proceed,  at  some  point  it  becomes  easier  to  form  a  new 
thrust  zont  within  the  shield  to  the  east  than  to  continue  uplifting  the 
mountains.  Suarez  et  al.  (1983J  envision  that  the  lateral  growth  of  the  Andes 
proceeds  by  successive  eastward  jumps  in  the  locus  of  thrusting  and  uplift. 
Such  outward  (eastward)  migration  of  foreland  thrusting  may  have  also  occurred 
in  western  North  America  [e.g.,  Armstrong  and  Oriel,  1965],  though  this  view 
is  not  held  by  all  [Mountjoy,  1966], 

Ue  suggest  a  scenario  for  the  growth  of  the  arc  in  the  Java-Bali  region  of 
the  Sunda  arc  that  is  similar  to  that  proposed  for  the  Andes  by  Suarez  et  al . 
[1983]  but  in  its  infancy.  The  Sunda  arc  example  also  suggests  that  discrete 
jumps  in  the  locus  of  faulting  proceed  along  strike  as  well  as  outward  from 
the  mountain  belt.  The  evolutionary  sequence  along  strike  is  exemplified  here 
by  the  spatial  and  temporal  development  of  the  Java,  Bali,  and  Flores  Basins. 

The  greater  widths  of  Java  and  Sumatra  than  those  of  the  islands  to  the  east 
are  not  due  to  voluminous  magmatic  activity  but  rather  because  of  the  exposure 
of  the  foreland  basins.  The  Java  Basin  in  the  Tertiary  was  a  shallow  marine 
basin  bordered  on  the  north  by  the  subaerial  Sunda  Shelf  and  to  the  south  by 
an  active  volcanic  arc  [Weeda,  1958],  (A  similar  basin  formed  in  Sumatra 
[Katili,  1974],  but  we  limit  our  discussion  to  the  Java  Basin.)  Prior  to  the 
development  of  the  Java  Basin,  the  Sunda  Shelf  probably  extended  to  the 
present  south  coast  of  Java  and,  accordingly,  the  Java  Basin  formed  as  a 
dovnwarp  in  the  Sunda  Shelf.  The  6-km-thick  sedimentary  sequence  in  the  Java 
Basin  reveals  a  marine  transgression  in  the  Oligocene,  bathyal  marine 
conditions  in  the  Miocene,  and  land  conditions  by  Pliocene  [Hamilton,  1 Q 7 9 ; 
Weeda,  1958],  Folding  in  the  west  Java  Basin  began  in  the  Miocene  when 
bathyal  conditions  prevailed  [Hamilton,  1979]  and  the  intensity  of  folding 
decreased  towards  the  Madura  basin  to  the  east  and  towards  the  volcanic  arc  to 


the  south  (Weeda,  1 9 5 8 J .  The  end  of  sedimentation  in  the  Java  Basin  is  marked 
by  emergence  and  a  major  orogenic  phase  in  middle  to  late  Pleistocene 
[Soetantri  et  al . ,  1973],  Today  in  the  Madura  and  Java  Basins,  the  younger 
sediments  are  less  deformed  than  the  older,  but  even  upper  Quaternary 
sediments  are  involved  in  the  compressional  deformation  [Hamilton,  1979; 

Weeda,  1958],  The  total  amount  and  present  rate  of  shortening  in  the  Java 
Basin  are  unknown  but  the  lack  of  earthquakes  or  gravity  anomalies  suggest 
that  they  are  small,  despite  the  evidence  for  continued  deformation  [Weeda, 
1958]  . 

Activity  began  in  the  Bali  Basin  as  it  waned  in  the  Java  Basin  and  we  are 
probably  witnessing  a  similar  progression  underway  in  the  Bali  Basin.  However, 
we  do  not  intend  to  infer  a  causal  relationship  between  the  onset  of 
deformation  in  one  basin  and  its  demise  in  the  other.  The  high  variability  in 
the  character  of  backarc  thrusting  at  present  along  the  Sunda  Arc  indicates 
that  the  thrusting  is  strongly  dependent  on  conditions  in  the  forearc  as  well 
as  the  backarc  [Silver  at  al.,  1983];  hence  the  eastward  migration  of  activity 
in  the  backarc  from  Java  to  Bali  may  reflect  some  forearc  process  (such  as 
migration  of  the  Roo  Rise),  as  easily  as  a  backarc  process  (for  instance,  the 
resistance  of  the  buoyant  Sunda  Shelf  to  thrusting  beneath  the  volcanic  arc). 
The  point  is  that  at  some  time  in  the  future,  the  Java  and  Bali  Basins  may 
appear  as  a  linear  zone  of  deformation  but  will  be  characterized  by  greatly 
different  histories.  Also  structurally  on  line  with  the  Bali  and  Java  Basins 
will  be  the  suture  along  which  the  Flores  Basin  closed.  The  Flores  Basin, 
while  undergoing  closure  concurrently  with  the  Bali  Basin,  is  morphologically 
distinct  from  the  Bali  Basin  in  that  it  is  floored  by  oceanic  crust  and  is 
much  deeper.  With  the  exceptions  of  subtle  sediraentological  facies  differences 
and  the  possibility  that  its  suture  will  be  marked  by  ophiolites,  the 
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structures  resulting  from  closure  of  the  Bali  and  Flores  Basins  may  be 
indistinguishable.  The  picture  will  be  complicated  further  if  new  fold  and 
thrusts  belts  develop  outboard  (north)  of  the  old,  as  in  the  scenario  of 
Suarez  et  al.  [1983]. 

Summary 

The  focal  depths  and  fault  plane  solutions  for  8  earthquakes  from  the  Bali 
Basin  have  been  determined  by  inversion  of  long-period  P  and  SH  waves  and 
short-period  P  waves.  The  results  indicate  that  thrusting  occurs  beneath  the 
southern  margin  of  the  Bali  Basin  at  depths  that  range  from  10  to  18  km  and  on 
planes  that  dip  to  the  south  at  angles  between  13’  to  35’ .  The  events  had 
seismic  moments  that  range  from  4  xlO17  to  7  xlO1*  Nm.  and  several  were  quite 
destructive.  Summation  of  seismic  moments  suggests  a  minimum  closure  rate  of 
4  ±2  mm/a  over  the  past  25  years. 

Characteristics  of  the  waveforms  and  the  earthquake  locations  with  respect 
to  crustal  models  constrained  by  gravity  and  seismic  data  suggest  that  they 
occurred  primarily  within  the  crystalline  basement,  beneath  the  sedimentary 
layer.  The  south  dipping  nodal  planes  that  we  infer  to  be  the  fault  planes  dip 
too  steeply  to  be  due  to  thrusting  on  a  zone  of  oecollement  between  the 
basement  and  the  sediment  layer  but  also  too  shallowly  to  be  analogous  to  the 
steep  faults  along  which  thick  basement  blocks  are  uplifted,  such  as  during 
the  Laramide  faulting  of  western  North  America  or  in  the  Parapeanas  Ranges  of 
the  Andes.  These  earthquakes  probably  reveal  some  intermediate  activity,  such 
as  basement  shortening  south  of  the  foreland  fold- and- thrust  belt  represented 
by  the  central  part  of  the  Bali  Basin. 
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The  similarity  of  the  slip  vectors  of  the  8  earthquakes  to  the  convergence 
direction  between  the  Asian  and  Indian  Ocean  plates  suggests  that  the 
thrusting  in  the  backarc  region  is  related  to  plate  convergence.  Collision  or 
strongly  coupled  subduction  is  indicated.  However,  because  the  arc  is 
extensively  cut  by  strike-slip  faults,  it  is  probably  quite  weak,  and  the 
thrusting  here  is  not  easily  relatable  to  the  collision  between  the  Australian 
continent  and  eastern  Sunda  Arc  that  occurs  several  hundred  kilometers  to  the 
east.  Moreover,  the  Mesozoic  age  of  the  subducting  lithosphere,  the 
morphology  of  the  forearc  south  of  Bali,  the  lack  of  thrust  events  beneath  the 
forearc,  and  a  preponderance  of  large,  shallow  normal -faulting  events  at  the 
Java  trench  suggest  that  the  upper  and  lower  plates  are  weakly  coupled.  The 
only  likely  scenario  involving  a  collision  is  interaction  with  the  Roo  Rise,  a 
broad  feature  with  approximately  1  km  of  relief  that  rises  above  the  Indian 
Ocean  floor  and  intersects  the  Java  Trench  south  of  Java  and  Bali. 

Finally,  we  propose  a  scenario  for  the  development  of  the  Bali,  Java,  and 
Flores  Basins  that  has  implications  for  the  development  of  continental  fold- 
and- thrust  belts  and  for  the  growth  of  continental  arcs  such  as  the  Andes  and 
western  North  America.  The  island  of  Java  has  nearly  doubled  in  width  during 
the  Tertiary  by  the  addition  of  a  wide  and  deep  sedimentary  basin  to  its 
backside.  This  basin  likely  formed  in  the  Sunda  Shelf  and  was  filled  by 
sediments  from  both  the  magmatic  arc  to  the  south  and  from  the  Sunda  Shelf  to 
the  north.  This  process  occurs  now  in  the  Bali  Basin.  Folding  and  thrusting 
within  the  Java  Basin  continue  to  the  present,  but  has  declined  in  intensity 
and  the  Java  Basin  is  now  exposed.  The  Bali  Basin  formed  east  of  the  Java 
Basin  and  is  now  in  the  process  of  being  closed  by  thrusting.  To  the  east  of 
the  Bali  Basin,  the  Flores  Basin,  of  oceanic  origin,  is  actively  closing  by 
thrusting  at  its  southern  margin.  This  example  suggests  that  individual  thrust 
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episodes  may  migrate  along  strike  as  well  as  outward  from  mountain  belts. 
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Appendix 


The  formal  uncertainties  in  the  earthquake  source  parameters  determined  by 
the  inversion  procedure  underestimate  the  true  uncertainties.  Here,  we  briefly 
assess  the  realistic  uncertainties,  examine  the  possible  trade-offs  between 
parameters,  and  present  our  strategies  for  deciding  on  the  best  solution.  The 
procedure  followed  is  to  fix  the  parameter  being  examined  at  values  that 
bracket  the  best  fitting  value,  and  solve  for  the  remaining  parameters.  In  the 
figures  that  follow,  only  selected  seismograms  are  shown  but  all  calculations 
were  done  with  all  the  available  seismograms  (Fig.  6). 

Focal  depths  For  the  range  of  depths,  source  durations,  and  mechanisms 
presented  in  this  paper,  events  2,  6,  and  8  are  probably  representative  and 
are  used  to  examine  depth  uncertainties.  For  event  2  (Fig.  Al),  the  strong 
trade-off  between  the  depth,  source  duration,  and  dip  angle  results  in  very 
little  increase  in  variance  between  5  and  13  km  depth.  At  depths  greater  than 
13  km  the  reflected  phases  are  too  late,  resulting  in  phase  mismatch  and  large 
variance.  We  feel  that  the  diagnostic  seismogram  is  at  AAE,  for  which  the 
shallower  solutions  give  a  poor  fit  to  the  P-wave  because  of  the  increase  in 
the  dip  angle.  From  this,  we  infer  that  event  2  is  no  deeper  than  16  km  and  no 
shallower  than  8  km.  Although  the  P-wave  seismogram  at  AAE  was  not  available 
for  event  1,  we  suggest  that  the  uncertainty  in  its  depth  is  also  ±5  km,  based 
on  the  similarity  of  the  seismograms  to  those  of  event  2,  and  the  proximity  in 
time  and  space  of  the  two  events. 

Event  6  (Fig.  A2)  shows  much  less  of  a  trade-off  between  the  depth  and  other 
parameters.  Although  the  source  duration  increases,  the  dip  angle  is  more 
stable  than  for  event  2  and  the  variance  increases  with  decreasing  depth.  This 
stability  in  the  dip  angle  is  due  to  the  increase  in  the  relative  strength  and 
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number  of  P-wave  seismograms  with  P  and  pP  phases  of  opposite  polarity  (such 
as  AAE  for  event  2).  For  the  shallower  assumed  sources,  the  initial  cycle  of 
the  calculated  seismogram  for  GUA  becomes  small  because  of  the  destructive 
interference  between  P  and  pP.  The  estimated  uncertainty  in  depth  of  ±4  km  for 
this  event  probably  also  applies  to  events  3,  4,  5,  and  7,  because  of  the 
similarity  in  their  mechanisms  and  station  distributions. 

Event  8  (Fig.  A3)  is  similar  in  depth  and  orientation  to  event  1  but  has  a 
longer  source  duration.  As  for  event  2,  there  is  a  trade-off  between  depth  and 
source  duration  but  over  a  greater  range  of  depths;  from  12  to  24  km.  Above  9 
km  the  variance  decreases  because  of  the  change  in  dip  but  again  this  produces 
an  inacceptable  fit  at  AAE. 

The  uncertainty  in  source  depths  and  the  variations  in  seismic  moment  with 
depth  introduces  an  uncertainty  of  approximately  50%  in  the  estimated  seismic 
moments  (Figs.  Al ,  A2 ,  and  A3).  An  additional  systematic  bias  in  the  depth 
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(not  included  in  Table  1)  is  caused  by  the  uncertainty  in  the  velocity  used  to 
calculate  seismograms.  If  the  average  velocities  in  the  earth  above  the  source 
range  from  5.5  to  7.5  km/s ,  the  bias  in  depth  will  be  about  ±20%  of  the 
calculated  depth. 

Orientation  The  uncertainties  in  strike,  dip,  and  rake  are  measured  with 
respect  to  the  N-dipping  nodal  planes.  For  event  2  (and  by  similarity,  for 
events  1  and  8),  the  uncertainty  in  strike  is  between  10  and  20*  (Fig.  A4) . 

For  this  mechanism  the  near  nodal  P  and  SH  waves  are  most  sensitive  to  strike 
but  there  is  a  trade-off  with  the  rake  angle.  The  uncertainty  in  strike  for 
event  3  (and  events  4-7)  is  greater  (±20-30*;  Fig.  A5)  than  that  for  event  2 
despite  the  decrease  in  trade-off  with  the  rake  angle. 

The  dip  angle  of  the  N-dipping  nodal  plane  is  constrained  within  ±5*  (Fig. 
A6).  The  most  obvious  constraints  are  the  polarities  of  the  P-waves  but  the  SH 
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wave  amplitudes  are  also  sensitive  to  the  dip  angle.  Because  the  dip  angle  is 
determined  largely  by  seismograms  at  the  northern  stations,  it  will  vary 
directly  with  the  take-off  angles  for  rays  to  these  stations.  For  a  range  in 
source  velocities  from  5.5  to  7.5  km/s,  the  take-off  angles  and  thus  dip 
angles  will  vary  by  ±5*  for  a  dip  angle  of  60*  and  by  ±2*  for  a  dip  angle  of 
75*  . 

The  uncertainty  in  the  rake  angle  for  event  3  is  approximately  20*  (Fig.  A7) 
and  is  probably  similar  for  the  other  events.  The  rake  angle  is  largely 
constrained  by  the  SH  waves  but  trades  off  strongly  with  the  strike  angle. 

Note  that  a  difference  in  the  rake  angle  of  20*  on  the  N-dipping  auxiliary 
plane  corresponds  to  differences  in  the  S- dipping  (fault)  plane  of 
approximately  30*  in  strike  and  less  than  5*  in  dip  (compare  in  Fig,  6  events 
3  and  5  whose  rake  angles  differ  by  18*). 

Throughout  the  tests  of  the  orientation  parameters,  the  depth  (±1  km), 
source  duration  (±1  s),  and  seismic  moment  (±10%)  remain  very  stable. 
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Figure  Captions 

1.  Tectonic  and  geographic  map  of  the  eastern  Sunda  Arc  and  vicinity.  Active 
volcanoes  are  represented  by  triangles  and  bathymetric  contours  are  in 
kilometers.  Thrust  faults  are  shown  with  teeth  on  the  upper  plate.  The 
dashed  box  encloses  the  study  area. 

2.  Bathymetric  and  free-air  gravity  maps  of  the  Bali  Basin  and  vicinity. 
Contours  are  in  meters  and  mGals  (10  ^  m/s2),  respectively,  and  are  dashed 
or  omitted  where  data  control  is  poor.  The  lighter  points  show  the  locations 
of  measurements.  Heavy  lines  show  the  locations  of  Lines  63,  67,  and  B22 ,  and 
the  brackets  enclose  the  data  used  to  construct  Line  1.  Large  dots  show  the 
locations  of  the  earthquakes  discussed  in  this  paper. 

3.  Free-air  gravity  profiles  and  line  drawings  of  seismic  reflection  profiles 
(from  Silver  et  al.  [1983])  across  the  Bali  Basin  near  Lombok.  The 


locations  of  these  lines  are  shown  in  Fig.  2. 

4.  Seismic  profile  crossing  the  southern  Sunda  Shelf  and  Bali  Basin  north  of 

Lombok  (from  Ben-Avraham  [1973]).  Location  of  profile  is  shown  in  Fig.  2. 

The  top  figure  shows  the  digitized  bathymetry  (dots)  at  high  vertical 

exaggeration  and  theoretical  curves  for  a  flexed  elastic  plate.  For  a  plate 

deflected  vertically  at  its  free  end,  the  solution  is  of  the  form 
•  Ax 

z(x)-z  cos(Ax)  e  ,  where  z  is  the  deflection  at  distance  x,  z  is  t.ie 
o  o 

deflection  at  x-0,  A-( (p^-p^) g/4D)k ,  g  is  the  gravitational  acceleration 

(9.8  m/s2),  and  D  is  the  flexural  rigidity  of  the  plate  [Watts  and  Talwani, 

1974].  Assumed  densities  for  the  mantle  and  w0.er  are:  p  -3300  kg/ra3 , 

m 

Pw-1000  kg/m3.  The  curves  are  for  flexural  rigidity  values  of  1,  5,  and  10 
xlO22  Nm  (as  marked). 
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5.  (a)  Seismicity  of  the  Bali  region  from  1964-1979  (all  events  listed  by  the 
ISC).  Depth  intervals  are  in  km.  (b)  Cross-section  of  seismicity  in  the  Bali 
region.  This  projection  passes  through  8*S,  115°E,  at  an  azimuth  of  10*.  and 
includes  all  events  within  300  1cm  and  recorded  by  more  that  50  seismograph 
stations.  In  both  plots,  triangles  show  the  locations  of  active  volcanoes. 

6.  Plots  of  waveforms  and  fault  plane  solutions  for  Bali  Basin  earthquakes. 

As  in  all  following  figures,  observed  seismograms  are  shown  by  solid  lines 
and  the  dashed  lines  are  calculated  seismograms.  On  the  left  in  each  plot 
are  the  P  waves  and  on  the  right  are  the  SH  waves.  Letters  at  the  lower 
right  of  station  codes  correspond  to  those  in  the  focal  sphere  and  an 
asterisk  indicates  that  the  seismogram  was  not  included  in  the  inversion. 

The  small  vertical  bars  show  the  time  window  used  in  the  inversion.  The 
normalized  source  time  function  is  shown  on  the  time  axis  and  the  amplitude 
scale  is  for  seismograms  (normalized  to  an  instrument  magnification  of  3000 
and  an  epicentral  distance  of  40*). 

7.  (a)  Lower  hemisphere  projections  of  fault  plane  solutions  near  Bali  Island 
determined  in  this  study.  Stippled  quadrants  are  those  with  compressional  P 
wave  first  motions.  P  and  T  axes  are  shown  by  large  dots.  Events  shown 
with  larger  focal  spheres  have  seismic  moments  greater  that  10ls  Nm. 

(b)  Cross-section  of  Bali  earthquakes  with  side  projections  of  their  fault 
plane  solutions.  The  dipping  lines  through  the  hypocenters  represent  the 
projections  of  the  inferred  fault  planes  and  the  vertical  lines  show  the 
range  of  possible  depths. 

8.  Short-period  P  wave  seismograms  for  events  3  and  4.  The  fault  plane 
solutions  and  source  structures  are  the  same  as  for  the  long-period  data 
(Fig.  6  and  Table  1).  Amplitudes  have  been  normalized  by  the  rms  amplitude 
within  the  inversion  window. 
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9.  Selected  short-  and  long-period  P  seismograms  for  event  3,  used  to  examine 
the  effects  of  the  assumed  source  structure  on  the  waveforms.  For  the  short - 
period  data,  the  amplitudes  are  normalized  (see  Fig.  8)  and  the  source 
orientation  is  fixed.  The  orientation  was  allowed  to  vary  for  the  long- 
period  data.  Clearly,  the  best  fit  to  both  the  long-  and  short-period  data 
is  for  case  E. 

10.  Structural  models  and  gravity  data  for  the  Bali  Basin.  The  numbers  in  each 
layer  give  its  density  in  kg/m3  (except  for  the  slab  in  which  the  density 
contrast  is  given),  (a)  Isostatically  balanced  model  along  Line  1.  The  top 
structural  model  is  for  the  crustal  structure  (vertical  exaggeration  of  4) 
and  the  bottom  shows  the  entire  assumed  structure  (no  vertical 
exaggeration),  (b)  Structural  model  beneath  Line  67  that  satisfies  gravity 
data  and  assumes  a  constant  thickness  for  the  crystalline  crust  of  the  Sunda 
Shelf.  The  small  lines  beneath  kms  250  to  280  show  the  relocated  positions 
of  the  earthquakes  (based  on  arrival  time  data)  relative  to  this  model.  The 
slopes  of  the  lines  show  roughly  the  dips  of  the  inferred  fault  planes  for 
the  events.  Note  that  the  waveform  data  suggest  that  the  earthquakes 
occurred  beneath  thick  sediments  that  are  probably  found  10  to  30  km  to  the 
north  of  the  positions  shown. 

Al .  Comparison  of  observed  (solid  lines)  and  calculated  (dashed  lines)  P  wave 
seismograms  of  event  2  for  a  range  of  depths.  The  depth  is  held  fixed  while 
the  remaining  parameters  (strike,  dip,  rake,  moment,  and  the  source  time 
function)  are  determined  by  the  inversion.  The  variance  is  normalized  by  the 
variance  of  the  best-fitting  solution.  Seismic  moment  is  in  units  of  10'® 

Nra . 

A2 .  Comparison  of  observed  and  calculated  P  wave  seismograms  of  event  6  for  a 
range  of  depths.  Format  as  in  Fig.  Al. 
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A3.  Comparison  of  observed  and  calculated  P  wave  seismograms  of  event  8  for  a 
range  of  depths.  Format  as  in  Fig.  Al. 

A4 .  Comparison  of  observed  and  calculated  P  and  SH  wave  seismograms  of  event  2 
for  a  range  of  strike  angles.  Format  as  in  Fig.  Al . 

A5 .  Comparison  of  observed  and  calculated  P  and  SH  wave  seismograms  of  event  3 
for  a  range  of  strike  angles.  Format  as  in  Fig.  Al. 

A6 .  Comparison  of  observed  and  calculated  P  and  SH  wave  seismograms  of  event  3 
for  a  range  of  dip  angles.  Format  as  in  Fig.  Al. 

A7 .  Comparison  of  observed  and  calculated  P  and  SH  wave  seismograms  of  event  3 
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Fig.  6.  Plots  of  waveforms  and  fault  plane  solutions  for  Bali  Basin  earthquakes  In 
all  plots,  observed  seismograms  are  shown  by  solid  lines,  and  the  dashed  lines  are 
calculated  seismograms.  On  Che  left  In  each  plot  are  the  P  waves,  and  on  the  right  are 
the  SH  waves.  Letters  at  the  lower  right  of  station  codes  correspond  to  those  in  the 
focal  sphere,  and  an  asterisk  Indicates  that  the  seismogram  was  not  Included  In  the 
Inversion.  The  small  vertical  bars  enclose  the  time  window  used  In  the  inversion.  The 
normalized  source  time  function  Is  shown  on  the  time  axis,  and  the  amplitude  scale  Is 
for  seismograms  (normalized  to  an  instrument  magnification  of  3000  and  an  epicentral 
distance  of  40*). 
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ABSTRACT 


An  overview  is  presented  of  the  stratospheric  cryogenic 
interferometer  balloon  experiment  (SCRIBE).  This  program  is 
designed  to  yield  high  resolution  emission  spectra  of  the 
earth's  atmosphere.  The  program's  history  is  briefly 
summarized.  The  data  encoding  and  transmission  schemes 
employed  are  described  and  the  procedures  used  for  data 
reduction  are  detailed.  The  results  of  data  analysis 
performed  by  groups  at  the  University  of  Denver,  the 
University  of  Massachusetts,  and  the  Air  Force  Geophysics 
Laboratory  are  reviewed.  A  bibliography  of  all  SCRIBE 
related  publications  and  tables  of  the  processed  spectra 
produced  by  the  instrument  are  included. 
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THE  SCRIBE  PROGRAM:  AN  OVERVIEW 


Introduction 

The  STRATOSPHERIC  CRYONGENIC  INTERFEROMETER  BALLOON 
EXPERIMENT  is  an  on-going  project  of  the  Optical  Physics 
Division  of  the  Air  Force  Geophysics  Laboratory  whose  goal 
is  to  obtain  high  resolution  atmospheric  emission  spectra  at 
up  to  40  km  altitudes.  These  measurements  are  of  importance 
to  the  Air  Force  because  the  emission  character isics  of  the 
stratosphere  must  be  considered  in  the  design  of 
surveillance  systems  which  look  at  targets  through  the 
stratosphere  or  against  a  stratospheric  background.  In 
addition  to  its  military  importance,  the  SCRIBE  data  may  be 
used  to  help  answer  a  variety  of  environmental  questions 
relating  to  the  role  anthropogenic  trace  gases  play  in  the 
photochemistry  of  the  atmosphere  and  their  effect  on  the 
earth's  radiation  budget. 

Initial  design  of  the  SCRIBE  instrument  began  in  the 
mid-1970's  with  flights  starting  in  1980  and  continuing 
until  the  present.  During  the  course  of  this  work  data 
reduction  and  analysis  have  been  performed  by  several 
groups,  using  similar  although  not  identical  methods.  Their 
results  have  not  been  collected  in  one  document  and  no 
updated  overview  of  the  SCRIBE  program  exists.  Because  of 
the  uniqueness  of  the  SCRIBE  measurements  and  the  increased 
interest  shown  in  them  by  researchers,  this  report  was 
prepared  to  provide  such  an  overview.  It  is  hoped  that  it 
will  serve  to  acquaint  potential  researchers  with  the 
program  and  the  data  which  it  has  generated  and  to  provide 
current  investigators  with  a  quick  reference  to  the  SCRIBE 
program . 
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The  SCRIBE  Instrument 


SCRIBE  is  a  balloon  launched  Fourier  tramsform 
spectrometer  ( FTS )  operated  at  liquid  nitrogen  temperature. 
Optical  alignment  over  a  wide  temperature  range  was 
accomplished  by  choosing  a  cat's  eye  interferometer  optical 
system  and  fabricating  all  mechanical  parts  from  the  same 
low  expansion  coefficient  material.  Mirror  position  is 
determined  by  a  mechanical  transducer  and  HeNe  laser 
fringes.  The  length  of  travel  of  the  moving  mirror  is  five 
cm.  The  FTS  is  mounted  in  a  large  dewar  equipped  wi'  an 
uncooled  ZeSe  window  of  low  emissivity.  The  entire  dewar  is 
tipped  to  allow  viewing  from  about  -5  degrees  to  +5  degrees 
about  horizontal.  A  flat  mirror  can  be  inserted  for  nadir 
viewing  and  an  ambient  temperature  blackbody  can  also  be 
inserted  to  aid  in  calibration.  A  gold  doped  germanium 
detector  operated  at  liquid  helium  temperature  is  employed. 
The  detector  signal  is  sampled  by  a  12-bit  analog  to 
digital  converter  whose  output,  along  with  a  variety  of 
housekeeping  signals  is  fed  to  a  custom  pulse  code 
modulation  (PCM)  encoder.  The  PCM  encoder  generates  a  72- 

4 

bit  frame  at  a  rate  of  10  frames  per  second.  The  data  is 
relayed  to  a  ground  station  via  an  S-band  telemetry  link. 
The  exact  form  of  the  PCM  stream  will  be  considered  in  more 
detail  in  the  discussion  on  data  reduction.  The  preceeding 
discussion  applies  to  the  instrument  as  it  was  flown  in 
1983,  modifications  for  a  specific  flight  are  described  in 
the  next  section. 


Program  Histor1 


In  1973  at  the  initiation  of  Dr.  George  Vanasse,  AFGL 
entered  into  a  contract  with  James  L.  Pritchard  of  Idealab, 
Inc.  for  the  development  of  a  cryogenical ly-cooled 
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interferometer.  At  the  conclusion  of  this  contract  in  1976, 
the  basic  SCRIBE  interferometer  as  it  exists  today  had  been 
constructed  [10,11].  In  its  initial  form  the  interferometer 
utilized  a  cooled  HgCdTe  detector  and  needed  major 
modifications  before  it  could  be  flown.  Under  a  contract 
from  AFGL  the  Department  of  Physics  of  the  University  of 
Denver  undertook  the  tasks  of  making  SCRIBE  flight  ready  and 
conducting  a  flight  program  [5].  Among  the  modifications 
they  made  were;  replacing  the  HgCdTe  detector  with  a  liquid 
helium  cooled  HgGe  detector;  interfacing  the  A/D  converter 
and  PCM  encoder;  repacking  all  electronics  to  ensure 
survival  at  40km. ;  mounting  the  dewar  in  a  gondola  with 
provisions  for  tilting  the  dewar  and  introducing  a  mirror  or 
blackbody  into  the  field  of  view.  On  October  8,  1980  SCRIBE 
was  flown  for  the  first  time  for  Holloman  Air  Force  Base  in 
New  Mexico.  Unfortunately,  problems  with  the  in¬ 
terferometer's  drive  and  loss  of  laser  fringes  limited  the 
data  obtained  to  a  few  noisy  spectrum  taken  at  low  altitude. 
The  instrument  was  flown  again  on  October  7,  1981  [17]  and 
June  16,  1982  but  in  each  case  problems  limited  the  value  of 
the  measurements  obtained  [6].  However,  these  three  flights 
served  to  demonstrate  the  potential  of  the  interferometer 
[13,14,20]  and  served  to  point  out  flaws  in  the  system  which 
were  subsequently  corrected  [12]. 

On  October  23,  1983  the  SCRIBE  instrument  as  described 
in  the  previous  section  was  successfully  flown  from  Holloman 
AFB ,  New  Mexico  [7].  All  systems  worked  well  during  the 
ascent  and  float  altitude  (  29  km}  was  reached  at  0750  MDT . 
At  some  point  the  balloon  developed  a  leak  and  after  30 
minutes  at  float  began  to  descend.  While  at  float,  a 
series  of  scans  at  various  look  angles  along  with  blackbody 
measurements  were  obtained.  During  the  descent  a  series  of 
nadir  viewing  scans  were  taken.  The  flight  was  terminated 
when  the  balloon  descended  to  18  km. 


Based  on  the  results  of  this  flight  an  optical  filter 
was  added  to  the  optical  path,  reducing  source  noise  from 
the  strong  15  micrometer  CO 2  band  and  resulting  in  improved 
performance  in  the  higher  wavenumber  region.  This 
instrument  was  flown  from  Roswell,  Texas  on  July  5,  1984  and 
reached  float  altitude  (30.8  km)  at  0745  MDT  [9].  All 
systems  worked  extremely  well  with  the  only  problem  being 
strong  winds  at  float  which  limited  the  time  the  balloon  was 
in  telemetry  range.  Good  data  was  obtained  on  ascent  and 
for  approximately  1  1/2  hours  at  float. 

After  the  1984  flight  the  primary  goal  of  obtaining 
high  resolution  atmospheric  emission  spectra  at  a  variety  of 
observation  angles  had  been  accomplished  and  it  was  now 
possible  to  plan  flights  designed  to  make  measurements  of 
specific  compounds  which  are  thought  to  play  key  roles  in 
atmospheric  photochemistry .  The  initial  compound  to  be 
studies  was  N^o^.  whose  presence  is  predicted  as  a  result  of 
reactions  between  NO  and  ozone.  The  interferometer  was 
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optimized  for  the  1250  cm  spectral  region  where  N2°5  hdS  a 
major  band  and  was  launched  on  June  21,  1985  at  0255  MDT  in 
order  to  be  at  float  just  before  dawn  when  the  N2°5  was 
predicted  to  peak  [9].  An  unexpected  wind  change  (just 
after  inflation  began)  made  for  a  difficult  launching.  Just 
after  release  the  payload  struck  the  ground  causing  a  loss 
of  command  and  balloon  control  links.  As  a  result, 
instrument  elevation  angle  and  detector  bias  could  not  be 
changed  and  the  data  from  this  flight  was  severely  degraded 
[4]. 

After  this  flight,  major  modifications  were  undertaken 
to  narrow  the  field  of  view  of  the  instrument,  provide 
better  control  of  the  viewing  angle,  and  visually  verify  the 
emission  source.  This  was  accomplished  by  adding  an 
uncooled  telescope  and  ground  controlled  steering  mirror  to 
the  optical  path  and  adding  a  bore-sighted  television  camera 


to  the  system.  An  operator  viewing  a  television  monitor 
could  control  the  view  by  manipulating  a  joystick.  The  new 
system  designated  SCRIBE-99  was  successfully  flown  on  August 
10,  1986.  The  balloon  was  launched  at  0700  MDT  and  reached 
a  float  altitude  of  76  Jcft.  All  systems  operated 
successfully  and  a  large  number  of  horizontal  and  nadir 
scans  were  acquired. 

Plans  for  future  flights  are  being  developed,  with  a 
high  priority  being  given  to  a  flight  plan  similar  to  the 
1985  flight. 

DATA  REDUCTION 

Because  the  large  volume  and  high  acquisition  rate  of 
data  by  the  SCRIBE  instrument  exceeds  the  capabilities  of 
flight  qualified  tape  recorders,  an  S-band  telemetry  link  to 
a  ground  recording  station  is  employed.  The  analog  output 
of  the  detector  preamplifier  is  fed  through  a  12-bit  analog 
to  digital  converter  whose  output  is  sampled  by  a  custom 
built  PCM  encoder.  An  interface  between  the  A/D  converter 
and  the  encoder  is  necessary  because  the  time  required 
between  inter ferogram  samples  is  slower  than  the  rate  at 
which  the  encoder  outputs  its  data  words.  Since  it  is 
crucial  that  a  sample  not  be  repeatedly  transmitted,  a 
signal  from  the  encoder  clears  the  A/D  outputs  and  sets  a 
data  valid  bit  to  zero.  Any  frames  sent  before  a  new 
measurement  is  made  contain  12  zero  data  bits.  When  a  new 
measurement  is  made  the  data  valid  bit  is  reset  to  one.  The 
12-bit  inter ferogram  sample  is  combined  with  four  status 
bits  to  form  a  16-bit  word.  Bit  13  is  tied  to  ground,  bit 
14  shows  drive  direction  (1;  forward),  bit  15  indicates  the 
detector  amplifier  gain  (1;  high  gain)  and  bit  16  indicates 
if  the  data  is  valid  (l;valid).  The  PCM  encoder  generates 
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72-bit  frames  at  a  rate  of  approximately  10^  frames/sec. 
Each  frame  consists  of  two  8-bit  sync  words,  an  8-bit 
counter,  the  16-bit  inter ferogram/ status  word,  two  8-bit 
words  indicating  carriage  position  and  speed,  and  two  8-oit 
words  containing  auxiliary  data.  Because  the  auxiliary 
measurements  do  not  require  as  high  a  sampling  rate  as  the 
interferometer,  they  are  commutated,  allowing  30  data 
channels  to  be  sampled  at  a  rate  of  700  samples  per  second. 
Among  the  measurements  on  these  channels  are  detector 
impendance,  amplifier  gain,  laser  fringe  amplitude,  various 
internal  and  external  temperatures,  elevation  angle, 
altitude  (pressure),  and  instrument  azimuthal  position. 
Other  data  was  recorded  on  these  channels  but  it  varied  from 
flight  to  flight.  The  University  of  Denver  group  has  the 
channel  assignments  as  well  as  the  calibration  data  which 
allows  these  digital  signals  to  be  interpreted.  On  the 
ground  the  telemetry  signal  after  passing  through  a  bit 
synchronizer  is  recorded,  along  with  millisecond  resolution 
time  code  (IRIGB)  data  on  analog  1/2"  magnetic  tapes  at  a 
speed  of  60  ips.  Simultaneously  the  output  of  a  parallel 
bit  synchronizer  is  directed  to  a  University  of  Denver  Nova 
1200  computer  acting  as  a  decommutator,  where  the 
inter ferogram  words  are  culled  from  the  data  stream, 
redundant  samples  discarded,  and  the  resulting 
inter ferograms  are  formatted  and  stored  on  digital  magnetic 
tapes.  Generally  at  least  two  ground  stations  make  analog 
recording  of  each  flight,  but  only  one  digital  recording  is 
done.  Ground  radar  stations  also  record  the  balloon's  X,  Y, 
Z  coordinates  and  radiosondes  are  flown  to  provide 
supporting  meteorological  data.  The  reduced  radar  and  sonde 
data  are  archived  at  AFGL . 

In  the  case  of  the  digital  tapes,  several  operations 
must  be  performed  before  they  can  be  transformed  into 
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spectra.  First  the  inter ferograms  must  be  corrected  for 
gain  changes  that  occur  during  and  between  scans.  These 
changes  are  necessary  because  of  the  widely  varying  radiance 
levels  observed  and  the  limited  resolution  of  the  12-bit  A/D 
converter.  The  first  gain  change  occurs  in  each 
inter ferogram  at  approximately  300  points  where  the  gain  is 
increased  by  a  factor  of  10.91  to  allow  the  region  beyond 
zero  path  difference  to  be  sampled  at  a  reasonable 
resolution.  During  the  flight/  the  peak  amplitude  of  each 
inter ferogram  is  monitored.  If  the  signal  falls  outside  the 
optimum  range  the  gain  is  automatically  adjusted  up  or  down 
by  a  factor  of  2  as  required  before  the  next  scan  is 
performed.  These  gain  changes  are  recorded  as  one  of  the 
auxiliary  data  channels  on  PCM  encoder.  Finally/  whenever 
the  instrument  views  the  blackbody  or  the  nadir/  the  gain  is 
reduced  by  changing  the  detector  bias  by  a  factor  of  2.1. 

After  the  gain  corrections  are  made,  a  spike  detector 
and  removal  program  is  run.  This  is  necessary  because 
occasionally  a  few  spurious  points  far  the  above  th?  normal 
signal  levels  are  seen.  They  are  believed  to  be  due  to  an 
electrical  problem  in  the  payload.  The  removal  program 
interpolates  new  values  at  these  points.  In  some  cases  a 
channel  spectrum  caused  by  the  window  is  present  and  must  be 
removed  in  a  similar  manner.  Next  a  program  is  run  which 
phase  corrects  the  inter ferogram.  A  fast  Fourier  transform 
is  done  to  produce  a  spectrum.  The  blackbody  spectra  can  be 
used  to  generate  an  instrument  function,  which  when  applied 
to  the  spectrum,  yields  a  calibrated  spectrum.  Using  this 
method  the  University  of  Denver  group  has  processed  selected 
spectra  from  the  1983  and  1984  flight.  From  the  1983  flight 
a  series  of  co-added  nadir  spectra  were  produced  (Table  1) 
and  from  the  1984  flight  four  sets  of  co-added  spectra  taken 


considerable  effort  this  was  traced  to  a  faculty  playback 
tape  recorder  which  was  repaired.  Presently  high  quality 
digital  tapes  can  be  produced  and  all  the  necessary  programs 
for  producing  finished  spectra  exist.  All  that  remains  is 
to  integrate  all  the  programs  to  allow  batch  processing  of 
the  large  volume  of  data  produced  by  the  SCRIBE  program. 

DATA  ANALYSIS 

Although  data  reduction  has  accounted  for  most  of  the 
effort  to  date,  data  analysis  has  resulted  in  several 
interesting  results. 

Dr.  Sakai's  group  found  features  aue  to  CO.,,  0g,  HNO  ^ , 
H2O  and  N2O  in  the  1983  spectra  [16]  The  Q  branch  of  the  l/^ 
band  of  CO2  at  667  cm"^'  which  is  very  opaque,  was  used  by 
them  to  measure  temperature  in  the  vicinity  of  the  balloon. 
The  values  obtained  were  about  10  K  higher  than  those 
measured  in  a  less  opaque  spectral  region  for  the  atmosphere 
at  a  similar  altitude  some  distance  from  the  instruments. 
They  attributed  this  discrepancy  to  a  mass  of  warm  air 
carried  along  with  the  balloon  on  its  ascent.  The  HNO 
features  were  identified  as  and  bands.  These  bands 

showed  maximum  emission  for  an  elevation  angle  with  a 
tangent  height  of  20  km,  indicating  that  the  concentration 
of  HNOg  is  a  maximum  around  the  tropopause.  By  a  similar 
analysis  it  was  concluded  that  the  majority  of  the  0g 
observed  occurred  below  30  km.  They  also  compared  that 
SCRIBE  spectra  to  synthetic  spectra  generated  from  the  AFGL 
line  computation. 

The  University  of  Denver  group  performed  a  similar 
analysis  [7]  and  in  addition  to  the  aforementioned  molecules 
also  indentified  features  due  CFCl^  (freon  11)  and  CF2CI2 
(freon  12).  They  also  obtained  a  laboratory  absorption 
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spectrum  of  HNO^.  By  comparing  this  spectrum  with  the 
SCRIBE  data,  the  assignment  of  H  NO  ^  features  in  the  SCRIBE 
spectra  was  verified.  In  addition,  the  laboratory  spectrum 
was  used  to  confirm  that  the  resolution  of  the  SCRIBE 
interferometer  is  0.06  cm  ^  { FWHM )  .  They  concluded  that  the 
instrument  has  sufficient  sensitivity  to  measure  atmospheric 
emission  at  high  altitudes  at  its  full  resolution.  Also  at 
the  Univerisity  of  Denver,  Dr.  Aaron  Goldman  prepared 
comparisons  of  simulated  and  SCRIBE  spectra  [3].  He  found 
generally  good  agreement  of  his  model  with  the  observed 
radiance,  but  noted  major  inadequacies  in  the  830  cm  ^ 
window  region  . 

Dr.  Goldman  has  also  prepared  similar  comparisons  for 
the  1984  flight  [1,2].  These  comparisons  were  made  for 
elevation  angles  of  1.9,  -3.2,  -3.7  and  90.0  degrees  while 
the  interferometer  was  at  float.  Again  reasonable  agreement 
was  found  but  a  number  of  adjustments  to  the  model  were 
necessary  before  this  could  be  accomplished. 

Dr.  Lawrence  Rothman  of  AFGL  has  compared  1983 
downlooking  spectra  obtained  at  61,  70,  81  and  88  kft  with 
synthetic  spectra  generated  by  the  FASCODE  modeling  program 
from  the  AFGL  line  parameter  file.  Preliminary  results 
indicate  that  there  is  a  major  discrepancy  between  the 
spectra  in  the  region  of  the  CO^  V  ^  Q  branch  at  667  cm  . 
This  discrepancy  is  not  yet  understood  and  additional 
analysis  is  being  undertaken  with  a  more  complete  set  of 
spectra . 

At  present,  a  more  detailed  analysis  of  the  1983  and 
1984  data  is  being  undertaken  by  personnel  from  the  AFGL  and 
OP ti Metrics,  Inc.  As  a  first  step,  an  atlas  of  the  spectra 
shown  in  Tables  2  and  5  is  being  prepared.  The  entire 


wavelength  region  of  each  spectrum  is  plotted  at  1  cm'1 
resolution,  and  the  opaque  region  around  667  cm"^  and  the 
window  region  around  810  cm  ^  are  plotted  at  0.06  cm'^ 


resolution.  Also  modifications  to  the  FASC0D2  computer 
L  model  are  being  done  to  allow  it  to  be  more  easily  used  in 
P  comparisons  with  the  SCRIBE  data. 

The  University  of  Denver  and  AFGL  are  also  beginning  to 


analyze  the  1986  data  but  no  results  are  currently 
available.  This  analysis  will  initially  focus  on  the  effect 
of  the  uncooled  telescope  on  quality  of  the  data  obtained. 
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is  referred  to  them  for  detained  descriptions  of  the  various 
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Table  1.  Co-Added  Spectra  From  the  October  23,  1983  SCRIBE 
Flight 


File  Co-Added  Altitude 

Name  Files  (kft) 


den . 

1 

Tape 

4 

Files 

86, 

87, 

88  , 

90,  91, 

92 

90.5 

den . 

2 

Tape 

5 

Files 

2, 

3,  4 

,  5, 

6,  7,  8 

,  9 

87.7 

den . 

3 

Tape 

5 

Files 

10, 

11, 

12, 

13,  14, 

15, 

16, 

17 

86.3 

den  . 

4 

Tape 

5 

Files 

18, 

19, 

20, 

21,  22, 

23, 

24, 

26 

85.5 

den . 

5 

Tape 

5 

Files 

27, 

28, 

29, 

30,  32  , 

33, 

35, 

36 

84.3 

den . 

6 

Tape 

5 

Files 

37, 

38, 

39, 

40,  41, 

42, 

43, 

44 

82.3 

den . 

7 

Tape 

5 

Files 

45, 

46, 

47, 

48,  49, 

50, 

51, 

52 

80.8 

den . 

8 

Tape 

5 

Files 

53, 

55, 

56, 

57,  58, 

61, 

62, 

63 

78.4 

den . 

9 

Tape 

5 

Files 

64, 

66, 

67, 

68  ,  69 , 

61, 

72, 

73 

75.0 

den . 

10 

Tape 

5 

Files 

74, 

75, 

76, 

77,  78, 

79, 

80, 

81 

70.0 

den . 

11 

Tape 

6 

Files 

2, 

3,  4 

,  5, 

6,  7,  8 

,  9 

61.0 

Notes  : 

1.  All  spectra  are  Nadir  views. 

2.  The  resolution  of  all  spectra  is  0.24  cm  1  FWHM. 

3.  A  blackbody  calibration  was  applied  to  all  spectra. 

4.  All  spectra  were  processed  by  the  University  of 
Denver  and  the  co-added  file  numbers  are  their 
designation . 

5.  Altitudes  are  approximate  because  the  balloon  was 
descending  during  this  period. 

6.  The  spectra  are  archived  on  magnetic  tape  and 
floppy  disk  at  the  AFGL . 


Table  2.  Calibrated  Spectra  Processed  by  the  University  of 
Massachusetts,  Amherst  from  the  October  23,  1983 
SCRIBE  Flight. 


File  Names 

Time 
(GMT  ) 

Altitude 

(kft) 

Elevation  Angle 
(degrees ) 

View 

SP830A1-A7 

12:16 

10 

7.5 

Horizontal 

M1-M7 

12:20 

20 

7.5 

Horizontal 

N1-N7 

12:25 

25 

7.5 

Horizontal 

01-07 

12:30 

30 

7.5 

Horizontal 

11-17 

13  :  15 

70 

1.7 

Horizontal 

21-27 

13:45 

95 

-90.0 

Nadir 

V1-V7 

13:50 

95 

-90.0 

Nadir 

W1-W7 

13:54 

95 

-90.0 

Nad  i  r 

B1-B7 

14:00 

95 

-90.0 

Nadir 

C1-C7 

14:05 

95 

-90.0 

Nadir 

D1-D7 

14:08 

95 

-90.0 

Nadir 

31-37 

14:13 

95 

-0.4 

Horizontal 

E1-E7 

14:17 

95 

-0.4 

Horizontal 

F1-F7 

14:20 

95 

-0.4 

Horizontal 

G1-G7 

14:24 

95 

*-0.4/2. 9 

Horizontal 

H1-H7 

14:27 

95 

-2.9 

Horizontal 

11-17 

14:  30 

95 

-2.9 

Hor izontal 

J1-J7 

14:33 

95 

*-2.9/-5.4 

Horizontal 

K1-K7 

14:37 

95 

-5.4 

Horizontal 

L1-L7 

14:40 

95 

-5.4 

Horizontal 

41-47 

14:45 

95 

-5.4 

Horizontal 

61-67 

14:48 

95 

-5.4 

Horizontal 

71-77 

14:53 

95 

- 

Black  Body 

81-87 

14:55 

95 

- 

Black  Body 

91-97 

15:02 

95 

7.5 

Horizontal 

P1-P7 

15:06 

95 

*7. 5/-90.0 

H/N 

Q1-Q7 

15:10 

95 

-90.0 

Nad  i  r 

R1-R7 

15:15 

90 

-90.0 

Nadir 

S1-S7 

15  :  18 

90 

-90.0 

Nadir 

T1-T7 

15  :  22 

90 

-90.0 

Nadi  r 

U1-U7 

15:25 

90 

-90.0 

Nadi  r 

Notes: 

1.  Resolution  of  all  scans  is  0.06  cm-  / FWHM . 

2.  *  Indicates  the  quantity  changes  during  the  interval. 

3.  A  Blackbody  calibration  was  applied  to  all  spectra. 

4.  Times  are  aprox imate  . 

5.  These  files  are  stored  on  magnetic  tape 

at  the  University  of  Massachusetts ,  Amherst. 


Table  3.  Tape  FTHY06 :  Calibrated  Spectra  from  the  October 
23,  1983  SCRIBE  Flight. 

File  Name  Time  Altitude  Elevation  Angle  View 
(GMT)  (kft)  (deg) 


SP830M4 

12:20 

20 

7.5 

Horizontal 

N 1 

12:25 

25 

7.5 

Horizontal 

01 

12:30 

30 

7.5 

Horizontal 

V3 

13:50 

95 

-90.0 

Nadir 

Wl 

13:54 

95 

-90.0 

Nadir 

B1-B7 

14:00 

95 

-90.0 

Nadir 

C1-C5 

14:05 

95 

-90.0 

Nadir 

Dl 

14:08 

95 

-90.0 

Nadir 

31,36 

14:13 

95 

-0.4 

Horizontal 

El 

14:17 

95 

-0.4 

Horizontal 

F6 

14:20 

95 

-0.4 

Horizontal 

G1-G3 

14:24 

95 

♦-0.4/-2.9 

Horizontal 

H1-H2 

14:  27 

95 

-2.9 

Horizontal 

11  ,  12 

14:30 

95 

-2.9 

Horizontal 

J2,  J7 

14:33 

95 

*-2. 9/-5 . 4 

Horizontal 

K3-K7 

14:37 

95 

-5.4 

Horizontal 

Ll ,  L2 

14:40 

95 

-5.4 

Horizontal 

41-47 

14:45 

95 

-5.4 

Horizontal 

95.97 

15:02 

95 

7.5 

Horizontal 

PI ,  P3 

15:06 

95 

*7.5/90.0 

Nadir 

01 

15:10 

95 

-90.0 

Nadir 

Notes  : 

1.  Resolution  of  all  scans  is  0.06  cm~^  / FWHM . 

2.  *  Indicates  the  quantity  changes  during  the  interval. 

3.  A  Blackbody  calibration  was  applied  to  all  spectra. 

4.  Times  are  aproximate. 

5.  These  files  were  prepared  by  Dr.  Sakai  at  the  University 
of  Massachusetts,  Amherst  and  are  stored  on  magnetic 
tape  at  AFGL 


Table  4.  Uncalibrated  Spectra  Processed  for  the  University 
of  Massachusetts,  Amherst  from  the  July  5,  1984 
SCRIBE  Flight 

Time 

(GMT)  File  Names  Elevation  Angle  Gain  Bias 

(deg  ) 


12:21-12: 22 

M I  -M  K 

4.9 

1 

2.1 

12:24-12:33 

QE , QF , QL , PE 
QJ  ,  QK  ,  P  F  ,  PC 

4.9 

2 

2.1 

12:50-13:45 

KA-KL , LA-LI  , 
LK, LL ,MA , MB , 
ME, HI ,HK , HL  , 
IA-IL , JA-JL 

1.9 

2 

2.1 

13:58-14:05 

GA , GB , G  D  —  G  J 

-0.6 

2 

2.1 

14:10-14:35 

H A-HH , TA-TH , 
UA-UH 

-3.2 

1 

2.1 

14:47-14:58 

AA-AL , BA , BB , 
BD-BF 

-3.7 

1 

2.1 

15:06-15:12 

DA-DH 

-1.2 

2 

2.1 

15: 28-15: 30 

EA-ED 

Blackbody 

1 

1 

15:45 

RB 

-1.2 

1 

1 

15:48 

RF 

-90.0 

2 

1 

15:50-15:53 

RH-RL 

-90.0 

1/2 

1 

Notes  : 

1.  Resolution  of  all  scans  is  0.06  cm  ~  . 

2.  No  blackbody  calibration  was  performed. 

3.  These  files  are  stored  on  magnetic  tape 

at  the  University  of  Massachusetts,  Amherst. 


Table  5.  Tape  FTHY85:  Uncalibrated  Spectra  from  the 
July  5,  1984  SCRIBE  Flight 


File  Name 

GMT 

Altitude 

(kft) 

Elev . Ang  . 

( deg ) 

Gain 

Bias 

847M  J 

12:22 

40 

4.9 

1 

2.1 

847QA 

12:45 

60 

1.9 

1 

2.1 

847  LC 

13:01 

7C 

1.9 

2 

2.1 

847HJ 

13:21 

75 

1.9 

2 

2.1 

847  I A 

13:24 

80 

1.9 

2 

2.1 

847  IK 

12:33 

85 

1.9 

2 

2.1 

847  J  B 

13:36 

90 

1.9 

2 

2.1 

847JI 

13:43 

94 

1.9 

2 

2.1 

847  FA 

13:49 

100 

1.9 

2 

2.1 

847GB 

13:58 

100 

-0.6 

2 

2.1 

847HB 

14:11 

100 

-3.2 

2 

2.1 

847  5C [TC ] 

14:22 

100 

-3.2 

2 

2.1 

847UC 

14:31 

100 

-3.2 

1 

2.1 

847AC 

14:48 

100 

-3.7 

1 

2.1 

847  BF 

14:58 

100 

-1.2 

1 

2.1 

846DC 

15:05 

100 

-1.2 

2 

2.1 

Notes  : 

1.  Resolution  of  all  scans  is  0.06  cm"^. 

2.  No  blackbody  calibration  was  applied  to  the 
spectra . 

3.  File  8475C [TC ]  was  obtained  from  a  coadd  of  eight 

interferograms. 

4.  These  files  were  prepared  by  Dr.  Sakai  at 
the  University  of  Massachusetts,  Amherst  and 
are  archived  on  magnetic  tape  at  AFGL . 
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Table  6.  Inter ferograms  Prepared  by  the  University  of 
Denver  from  the  July  5,  1984 


Tape 

Files 

Elevation  Angle 

Altitude 

Gain  Bias 

Number 

( deg ) 

(kft) 

1 

2-66 

4.9 

* 

1/2 

1 

1 

67-100 

4.9 

★ 

1 

2.1 

1 

101-125 

4.9 

★ 

2 

2.1 

2 

2-6 

4.9 

★ 

2 

2.1 

2 

12-113 

1.9 

* 

n 

2.1 

3  A 

2-74 

1.9 

★ 

2 

2.1 

3  A 

78-112 

-0.6 

100.3 

2 

2.1 

4A 

2-10 

-3.2 

100.3 

2 

2.1 

4  A 

11-73 

-3.2 

100.3 

1 

2.1 

4A 

79-98 

-3.7 

100.3 

1 

2.1 

5  A 

2-39 

-3.7 

100.3 

1 

2.1 

5  A 

43-97 

-1.2 

100.3 

2 

2.1 

6A 

5-7 

-1.2 

100.3 

2 

2.1 

6A 

11-46 

Blackbody 

100.3 

1 

1 

6A 

51-83 

-90.0 

100.3 

1/2 

1 

6A 

87-96 

-1.2 

* 

* 

* 

Notes : 

1.  *  indicates  the  quantity  varied  during 
the  interval. 

2.  These  inter ferograms  are  stored  on  magnetic 
tape  at  AFGL . 

3.  Files  are  not  listed  where  the  data  was 
questionable,  for  example,  when  the 
interferometer  was  being  moved  to  a  new 
elevation  angle. 
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|  ABSTRACT 

I 

■  Modifications  to  the  previously  created  one-dimensional,  periodic  plasma 

;  simulation  code  are  described.  These  modifications  allow  simulation  of 

diode-type  systems  in  which  the  particles  are  not  periodic  and  various 
boundary  conditions  nay  be  imposed  on  the  fields.  Particles  are 
\  injected  from  reservoirs  on  either  side  of  the  simulation  box  according 

;  to  some  predetermined  maxwellian  distribution  function;  particles  that 

I  exit  the  system  on  either  side  are  absorbed  by  that  wall.  The  walls 

!  function  as  capacitor  plates:  when  a  particle  leaves  the  wall  it  leaves 

behind  a  charge  equal  and  opposite  to  its  own,  and  when  the  particle  is 

|  absorbed  by  a  wall  it  contributes  its  charge  to  that  wall.  The 

]  potential  in  the  box  is  determined  by  Poisson's  equation  subject  to  one 

i 

|  of  three  possible  boundary  conditions:  1)  the  potential  is  set  equal  to 

zero  at  both  ends,  2)  the  potential  is  fixed  at  zero  on  the  left  and  a 
specified  value  on  the  right,  or  3)  the  potential  is  zero  on  the  left 
1  and  floats  on  the  right.  In  case  3  the  charge  in  the  walls  makes  a 

difference.  Tests  of  the  nonperiodic  part  of  the  simulation  code  are 
described. 
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I.  INTRODUCTION 


A  plasma  simulation  with  periodic  boundary  conditions  on  the  fields 
and  particles  simulates  infinite  plane  waves  in  an  infinite  medium.  If 
one  starts  one  species  of  particles  with  a  given  drift  and  then  lets  the 
system  develop,  one  is  watching  the  decay  of  an  infinite  beam  as  it 
propagates  through  the  infinite  medium.  If  one  is  instead  interested  in 
the  evolution  of  plasma  in  one  spatial  location  through  which  a  beam  is 
directed,  a  nonperiodic  system  is  more  appropriate.  A  nonperiodic 
system  is  also  useful  in  modeling  various  configurations  of  isolated 
probes  in  a  surrounding  plasma.  Since  one  of  the  principal  interests  at 
this  research  location  is  the  propagation  of  beams  through  plasma,  the 
development  of  the  simulation  code  was  continued  so  as  to  allow 
simulation  of  nonperiodic  systems. 

The  simulation  code,  except  for  the  nonperiodic  options,  is  a  subset 
of  ESI  (see  Birdsall  and  Langdon,  1985).  The  momentum  conserving, 
linear  weighting  code  is  used  without  any  provision  for  smoothing  of 
charge  densities  or  fields.  As  in  ESI,  field  variables  are  calculated 
on  a  grid  and  linearly  interpolated  to  particle  positions.  The 
particles  themselves  are,  in  effect,  finite  width  planar  charge  sheets. 
In  the  periodic  part  of  the  code,  periodic  boundary  conditions  are 
imposed  on  both  the  fields  and  the  particles;  when  a  particle  exits  one 
end  of  the  simulation  box  it  is  reinserted  in  the  other  end  with  the 
same  speed.  One  may  study  propagation  of  waves  at  anv  angle  with 
respect  to  a  uniform  magnetic  field.  Diagnostm  data  sm  h  as  eiectri' 
field  profiles,  phase  space  information,  and  information  ui  part  r r le 
distribution  functions  ( much  the  same  as  m  FS l '  are  wr 1 1 1  en  t  i 


separate  data  file  which  can  bn  used  t,v  sepirite  r*.  ’  <  pi  *?  t  ‘m 


diagnostics.  One  of  these  programs  just  plots  the  data  in  the  order  it 
was  written.  The  other,  a  postprocessor,  can  plot  any  interval  of  time 
of  history  data,  and  can  also  make  and  plot  discrete  Fourier  transforms 
of  these  time  intervals. 

Since  the  possibility  of  nonperiodic  simulations  was  always  in  mind, 
provisions  were  made  in  the  original  coding  to  aid  in  such  simulations. 
In  this  latest  research  period,  these  provisions  were  used,  modified, 
and  expanded.  Modifications  made  to  allow  for  nonperiodic  particles, 
and  modifications  in  calculating  the  fields  in  nonperiodic  systems  are 
described  below. 

II.  OBJECTIVES 

The  objectives  of  this  research  effort  were: 

1)  To  continue  development  of  the  nne-dimensir  nal ,  electrostatic 
simulation  code, 

2 1  To  update  the  user  document  at ;or  and  educate  people  in  the  use  of  the 

c  ode , 

<i  t  ■,  us,»  the  '  de  : sup:  -r*  t  lata  analysis  *  r  >m  beam-plasma 
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of  this  material  are  given  in  the  users  guide  which  has  been  prepared 
for  workers  on  site. 

System:  The  simulation  system  may  be  thought  of  as  a  diode  with 
capacitor  plates  for  its  two  walls.  One  may  think  of  the  plates  also  as 
particle  guns  that  leak  or  direct  plasma,  pulled  from  the  plates 
themselves,  into  the  system.  The  two  plates  may  either  be  joined  bv  a 
conductor  (short  circuit  conditions),  connected  to  a  battery  (fixed 
potential  drop),  or  left  to  float  (no  connection).  The  picture  is 
complicated  somewhat  by  the  option  of  treating  any  given  species  as 
periodic.  In  this  case,  the  periodic  species  sees  the  system  as  a 
toroid;  when  the  particle  exits  one  end,  it  enters  the  other  end  with 
its  velocity  unchanged. 

Particles:  To  accommodate  the  injection  of  particles,  provisions 
were  made  for  creation  of  reservoirs  from  which  particles  are  injected 
at  each  time  step.  The  reservoirs  are  loaded  with  Maxwellian 
distributions  in  the  same  manner  that  particles  were  loaded  in  the 
periodic  code,  and  thereafter  are  not  changed;  one  cycles  repeatedly 
through  the  list  of  reservoir  velocities  as  needed  to  choose  the 
velocities  of  entering  particles.  In  the  reservoirs,  particles  are 
sorted  into  left  going  and  right  going  particles,  and  the  initial 
numbers  of  particles  in  each  of  these  groups  are  recorded.  One  has  the 
option  of  setting  a  cutoff  number  such  that  if  the  number  of  right  or 
left  going  particles  falls  below  that  number,  no  particles  of  that  group 
>•*'  injected. 

•'cuiib  particles  are  injected  each  timestep  to  reflect  the  proper 

,  •  «>n '  <•  r  :  ug  particles  implied  by  the  Maxwellian  distribution  and 

"  ■•  ••■  -‘‘isen  fur  that  particle  species.  In  practice,  the  flux  is 
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Fields:  The  first  change  from  the  periodic  code  is  ’n  how  the  charge 

density  is  calculated.  In  both  the  periodic  and  nonperiodic  codes,  the 
charge  due  to  a  particle  is  linearly  distributed  between  the  two 
neighboring  grid  points  around  the  particle.  In  the  periodic  code,  the 
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*.  r  •  !  pom's,  ’  he  '  r  :  d  l  agona !  solver  dps<  r  l  f>p  1  in  appendix  I1  if 
i'.rfsall’s  bo  oe  hirfsall  and  l.angdon  i'-fH'i)  is  employed.  If  float  mu 

r  i  a’:*  hard  boundary  is  assumed,  *  he  potent  t.il  is  tabulated  via  a 
en'»rei!  difference  f  ;rti  of  Poisson's  equation  with  the  boundary 
■'ill*  is  r  rhat  zero  'barge  resi  les  out  side  grid  points  1  and  NO  *■  1,  and 

•ha'  the  potential  of  grid  point  I  is  zero.  The  calculation  proceeds 

f  rotn  left  to  r  i st h r  .  If  short  <  in mt  boundary  conditions  are  assumed,  a 
linear  term  is  added  to  the  potential  at  each  grid  point  found  from  the 
’floating*  calculation  so  that  the  potential  of  grid  point  NG  +  1  is 
zero.  In  all  cases,  the  electric  field  on  the  two  boundaries  is 
calculated  via  a  two  point  difference  formula,  and  the  electric  field 
within  is  calculated  via  a  centered  difference  formula. 

Tests :  Except  for  short  runs  to  debug  the  code,  only  two  types  of 
simulations  have  been  run  which  might  serve  as  tests  of  the  nonperiodic 
code.  These  are  1)  an  attempt  to  reproduce  simulations  run  in  Chodura 
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(1982)  involving  injection  of  unequal  fluxes  of  electrons  and  ions 
through  one  boundary  into  an  absorbing  wall  at  the  other  boundary  (with 
no  magnetic  field  present),  and  2)  a  two-stream  simulation  in  which  cold 
electrons  of  equal  density  are  injected  from  opposite  walls.  The  intent 
here  is  to  compare  the  results  of  the  nonperiodic  system  to  those  of  the 
periodic  simulation  of  a  similar  system. 

The  results  are  inconclusive,  but,  in  case  2,  promising.  In 
Chodura's  simulations  a  steady  state  was  arrived  at  in  which  the 
potential  made  a  smooth  transition  from  zero  at  the  injection  boundary 
to  a  negative  potential  at  the  absorbing  wall.  The  wall  potential 
depended  on  a  straightforward  way  on  the  ratio  of  the  injected  fluxes. 

In  simulations  run  with  the  present  code,  no  steady  state  was  arrived  at 
even  after  2000  plasma  periods.  Even  potential  profiles  averaged  over  4 
plasma  periods  (  251  time  steps)  showed  no  sign  of  settling  down.  The 
wall  potentials  in  the  averaged  profiles  were  consistently  below 
(sometimes  below  half)  as  much  as  the  theory  predicted.  It  could  be  the 
injection  scheme  that  is  at  fault  (too  much  noise?),  or  something  else. 

In  the  case  of  the  two-stream  simulation,  the  results  seem  more 
straightforward.  Figures  1-5  show  a  comparison  between  diagnostics 
taken  in  the  periodic  case  (left  column  in  each  Figure)  and  those  taken 
in  the  nonperiodic  case  (right  column  in  each  Figure).  In  both  runs, 
the  simulation  boxes  were  initially  filled  with  128  particles  of  each 
species  with  the  following  parameters: 

L-2fT,  DT-0.2,  NG-32  for  the  system,  and  N-128,  UT-l.O,  QM— 1 .  ,  VO-  1, 
and  an  initial  perturbation  of  position  according  to 
x  — >  x  +  0.01cos( 2*x) .  Similar  behavior  is  observed  in  both 
simulations,  with  differences  that  one  might  expect.  For  instance,  the 
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mode  that  was  initially  excited  (mode  2)  grows  in  both  cases  with  about 
the  same  growth  rate  (see  Fig.  1).  The  phase  space  of  right  going 
particles  (Fig.  2)  and  right  going  particles  (Fig.  3)  have  similar 
features,  and  the  growth  of  mode  2  is  clearly  seen  initially.  Also  the 
potential  and  electric  fields  (Figures  4  and  5)  show  similar  features. 
However,  as  can  be  seen  in  the  phase  space  diagrams,  the  symmetry  of  the 
periodic  phase  space  pictures  between  species  one  and  two  (right  and 
left  going),  and  for  each  species  individually,  is  lacking  in  the 
nonperiodic  diagrams.  This  is  to  be  expected  since  the  velocity  of  the 
nonperiodic  particles  is  clamped  to  plus  or  minus  one  at  one  boundary. 
The  potential  of  the  nonperiodic  simulation  eventually  deviates  from 
that  of  the  periodic  case  because  of  reflection  of  particles  in  the  box, 
and  absorption  of  particles  by  the  walls;  note  that  the  integral  of  the 
potential  with  respect  to  position  is  not  zero  in  the  final  panel  The 
general  correspondence  between  these  two  simulations  in  the  early 
development,  however,  leads  to  the  qualitative  conclusion  that  there  is 
good  physics  represented  here.  Without  a  detailed  model  of  what  one 
should  expect  in  the  nonperiodic  case,  however,  one  cannot  make  more 
quantitative  statements. 

IV.  Documentation  of  Code 

The  users  guide  has  been  updated  to  reflect  the  new  changes,  and  to 
describe  how  to  add  diagnostics  to  the  code  at  will.  Many  comments  have 
been  added  to  the  source  listings.  In  addition,  several  people  have 
been  instructed  in  the  use  of  the  code. 
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V.  Support  of  Data  Analysis 


Only  one  simulation  with  inconclusive  results  was  run  to  simulate  a 
beam  passing  through  a  plasma.  The  bulk  of  the  time  was  occupied  with 
objectives  1  and  2. 

IV.  Recommendations 

a)  Given  the  users  guide,  a  copy  of  the  source  programs  on  a  floppy  disk 
and  on  tape,  and  the  instruction  of  a  few  people  in  the  use  of  the  cole, 
the  simulation  code  should  be  readily  available  and  useful  after  I  have 
left.  When  the  operating  system  of  the  lo  al  computer  changes,  some 
changes  may  have  to  be  made  in  the  plotting  routines;  one  may  wish  to  <i  > 
so  any  way  in  order  to  get  better  graphs.  The  programs  should,  however, 
be  useful  as  they  stand. 

b)  Better  tests  should  be  devised  for  the  existing  code.  There  should 
be  a  wealth  of  literature  on  the  Pierce  diode  (see  Birdsall  book  for 
reference  or  the  PDW1  reference  guide)  from  which  one  might  extract 
useful  test  simulations.  Another  possible  test  would  be  to  calculate 
the  linear  evolution  of  the  two  stream  system  pictured  in  the  Figures 
here,  and  to  compare  that  evolution  with  what  was  seen. 

c)  One  modification  in  the  particle  injection  scheme  may  be  worth 
trying.  Instead  of  randomly  locating  the  incoming  particles,  a  smoother 
injection  scheme  would  locate  the  particles  at  injection  sites  spaced  a 
distance  L/N  apart  and  travel  from  the  boundary  at  an  average  velocity 
determined  by  the  specified  distribution  function.  In  the  case  of 
injecting  cold  particles,  sites  at  intervals  of  1,/N  apart  would  flow 
into  the  simulation  box  at  a  speed  equal  to  that  of  the  initial  drift. 
The  added  complexity  might  be  repaid  with  a  quieter  simulation. 
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